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Part  I 


Introduction 


PART  I  -  INTRODUCTION 


1.0  General 

Intelligent  material  systems  and  structures  have  been  investigated  by  various  government 
agencies  and  private  corporations  in  the  past  decade.  Intelligent  material  systems  and 
structures  are  engineered  systems  with  a  significant  degree  of  autonomy  that  enables 
them  to  maximize  performance,  provide  adaptive  functionality,  and  minimize  life-cycle 
cost.  These  engineered  systems  are  formed  by  incorporating  highly  integrated  actuators 
and  sensors  into  the  structure.  These  components  may  have  structural  functionality,  as 
well  as  control  logic,  signal  conditioning  and  power  amplification  electronics.  Such 
actuating,  sensing,  and  signal  processing  elements  are  incorporated  into  a  structure  for  the 
purpose  of  influencing  its  mechanical  and  thermal  characteristics  to  achieve  desired 
structural  performance 

A  wide  variety  of  applications  for  intelligent  structure  technologies  exist  in  aerospace, 
defense,  auto,  civil,  and  medical  devices  industries.  These  include  aeroelastic,  control  and 
maneuver  enhancement,  reduction  of  vibrations  and  structure  borne  noise,  jitter  reduction 
in  precision  pointing  systems,  shape  control  of  plates,  trusses  and  lifting  surfaces, 
isolation  of  offending  machinery  and  sensitive  instruments,  flexible  robotic  arms  damage 
detection,  reduction  of  life-cycle-cost  of  buildings,  bridges  and  lifelines,  and  sensitive 
biocompatible  medical  instrumentation. 

The  rebuilding  and  enhancement  of  our  Nation’s  infrastructure  can  greatly  profit  fi-om  the 
research,  development,  and  implementation  of  intelligent  material  systems.  Intelligent 
material  systems  can  allow  us  to  design  buildings,  bridges,  lifelines,  etc.  to  not  only  have 
a  greater  envelope  of  utility  but  also  do  so  while  reducing  the  overall  life-cycle  cost — it  is 
precisely  this  objective  which  defines  intelligent  systems.  Intelligent  structures 
incorporate  irmovative  material  compositions  and  sophistication  in  the  design  and 
architecture  in  order  to  simplify  construction,  reduce  maintenance,  and  decrease  cost  of 
the  entire  system.  The  synergy  of  actuators,  sensors,  and  controls,  will  not  only  allow 
structures  of  the  future  to  have  increased  performance  and  functionality,  but  will  facilitate 
the  solution  of  many  socio-economic  issues  concerning  the  repair,  enhancement  and 
construction  of  our  Nation’s  infi'astructure. 

Several  technical  developments  have  been  combined  to  establish  the  potential  feasibility 
of  intelligent  structures.  The  first  is  a  transition  to  laminated  materials.  In  the  past, 
structures  were  built  fi'om  large  pieces  of  monolithic  materials  which  were  machined, 
forged,  or  formed  to  a  final  structural  shape,  making  it  difficult  to  incorporate  any  active 
element  inside  the  structure.  However,  in  the  past  thirty  years  a  transition  to  laminated 
material  technology  has  occiured.  Laminated  materials  allow  for  the  easy  incorporation 
of  active  elements  within  the  build  up  of  structural  forms. 
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Exploitation  of  the  off-diagonal  terms  in  the  material  constitutive  relations  is  a  second 
trend  makes  the  development  of  intelligent  structures  promising.  The  full  constitutive 
relations  of  a  material  include  characterizations  of  its  mechanical,  optical, 
electromagnetic,  chemical,  physical,  and  thermal  properties.  Researchers  have  focused 
on  only  block  diagonal  terms  for  quite  sometime.  For  example,  those  interested  in 
exploiting  a  material  for  its  structural  benefits  have  focused  only  on  the  mechanical 
characterization,  and  those  interested  in  exploiting  its  electrical  properties  have  focused 
on  the  electrical  characterization.  However,  much  gain  can  be  obtained  by  exploiting  the 
off-block  diagonal  terms  in  the  constitutive  relations,  which  for  example,  couple  the 
mechanical  and  electrical  properties  or  couple  the  mechanical  and  thermal  properties. 
The  characterization  and  exploitation  of  these  off-diagonal  material  relations  had  led  to 
many  of  the  progress  in  the  intelligent  structure  technology. 

Central  to  the  emergence  of  intelligent  structures  is  the  development  of  information 
processing,  artificial  intelligence,  and  control  disciplines.  These  developing  technologies 
have  created  the  enabling  infi'astructure  based  on  which  intelligent  structures  can 
develop.  Control  architectures  and  algorithms  are  the  brains  of  intelligent  systems.  These 
are  where  the  very  intelligence  comes  from.  The  design  of  these  architectures  and 
algorithms  is  of  utmost  importance  for  the  effective  utilization  of  smart  functions  of 
intelligent  structures. 

In  an  intelligent  structure,  the  presence  of  active  elements  such  as  actuators,  sensors  and 
processors  impact  the  host  structure  by  increasing  the  mass  and  stiffiiess  of  the  system 
and  interfering  with  the  load  path  and  potentially  introducing  new  structural 
discontinuities  which  must  be  accommodated.  This  introduces  great  challenges  for 
accurate  modeling  of  these  systems.  On  the  other  hand,  the  complex  interactions  among 
sensors,  actuators,  the  host  materials,  and  the  processors  must  be  fully  characterized 
before  the  technology  can  reach  its  full  potential.  Experimental  models  are  powerful  tools 
for  characterizing  these  complicated  interactions  and  showing  the  dynamic  behavior  of 
laboratory  scale  active  systems.  However,  experimental  models  are  limited  by  size,  cost, 
time,  noise,  and  many  laboratory  unknowns.  It  is  virtually  impossible  to  generate  large 
quantities  of  data  for  full  size  structures.  Therefore,  development  of  a  simulation 
package  for  design  and  analysis  of  smart  structural  systems  would  greatly  enhance  this 
technology. 

The  present  research  is  directed  toward  modeling  and  simulation  of  some  of  the  elements 
of  smart  structures  such  as  sensors,  actuators  and  host  structures.  Shape  Memory  Alloys 
(SMA)  are  among  candidate  materials  suitable  for  wide  applications  as  sensors  and 
actuators  in  smart  structures.  They  have  the  unusual  material  property  of  being  able  to 
sustain  and  recover  large  strains  (of  the  order  of  10%)  without  inducing  irreversible 
plastic  deformation  and  to  “remember”  a  previous  configuration  and  return  to  it  with  a 
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temperature  change.  Their  ability  to  achieve  large  strains  instantaneously  enables  the 
design  of  structures  capable  of  extremely  large,  recoverable  deflections.  In  addition, 
shape  memory  alloys  are  relatively  lightweight,  biocompatible,  easy  to  manufacture  and 
have  a  high  force  to  weight  ratio  [Wayman,  1980].  These  characteristics  of  SMA  make  it 
an  ideal  candidate  for  shape  control  and  frequency  tuning  of  composite  laminates.  On 
the  other  hand,  laminated  composites  are  best  suitable  as  host  material,  for  smart 
structural  systems.  Their  laminated  structures  provide  the  capability  of  embedding 
sensors,  actuators,  and  the  electronic  circuitary  in  the  load-bearing  structure.  Following 
is  a  summary  of  activities  related  to  the  technology  of  smart  structures  that  was  pursued 
during  the  course  of  this  research. 


2.0  Shape  Memory  Alloys  (SMA) 

Given  the  variety  of  potential  uses  for  SMA  and  the  high  interest  in  developing  new 
applications,  the  ability  to  accurately  model  and  analyze  structures  containing  SMA 
components  via  a  finite  element  procedure  is  extremely  attractive.  The  incorporation  of 
the  SMA  finite  element  procedure  into  the  design  stages  of  new  products  could  reduce 
development  times  and  costs  dramatically.  Since  the  properties  of  a  particular  alloy  can 
be  easily  and  drastically  altered  in  the  manufacturing  process,  the  properties  of  a  SMA 
component  for  a  given  design  can  be  varied  systematically  in  the  finite  element  analysis 
before  production.  This  optimization  procedure  will  enable  use  of  shape  memory  alloy 
components  with  specifically  tailored  properties  that  will  realize  their  full  potential  in 
each  individual  application. 

The  following  research  work  related  to  modeling  of  Shape  Memory  Alloys  (SMA)  were 
conducted  during  this  project. 

•  Development  of  thermomechanical  constitutive  relations  for  SMA  based  on 
the  kinetics  of  solid-solid  phase  transformation  for  arbitrary  paths  in  stress- 
temperature  space. 

•  Development  of  criteria  for  complete  cyclic  loading-imloading  and  partial 
loading-unloading  conditions. 

•  Development  of  numerical  schemes  in  a  finite  element  framework  for  state 
determination  during  phase  transformation  and  solution  of  highly  nonlinear 
field  equations  for  SMA. 

•  Implementation  of  nonlinear  truss  and  beam  elements  for  SMA. 

•  Introduction  of  a  3D  model  based  on  the  extension  of  the  proposed  ID 
thermomechanical  model. 

•  Formulation  of  large  kinematics  for  pseudoelastic  behavior  of  SMA  and  its 
finite  element  implementation. 
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•  Extension  of  the  thermomechanical  model  to  tension-compression  through 
introduction  of  new  internal  variables. 

•  Numerical  and  algorithmic  developments  with  respect  to  the  finite  element 
method  for  a  multi-dimensional  truss  and  a  two-dimensional  beam  elements. 

•  Extension  of  above  formulation  to  finite  kinematics  elements. 

Computational  aspects  of  solid-solid  phase  transformation  modeling  were  a  challenge 
during  the  course  of  this  research.  State  determination  with  all  possibilities  in  a  stress- 
temperature  space  required  design  of  specific  numerical  schemes  for  convergence.  The 
hysteresis  of  phase  transformation  in  cyclic  loadings,  its  expansion  and  contraction 
during  phase  transformation,  and  the  nonlinear  material  response  within  the  hysteresis 
loop  made  the  modeling  phase  transformation  in  Shape  Memory  Alloys  much  more 
complicated  as  compared  to  the  modeling  of  other  material  nonlinearities  such  as 
plasticity. 

The  details  of  these  developments  are  discussed  in  Part  II  of  the  present  report  and  in  the 
progress  reports  during  the  course  of  this  research. 


3.0  Multilayer  Shell  Theory  of  Laminated  Composites 

In  the  last  two  decades,  composites  have  found  increasing  application  in  many 
engineering  structures.  Recent  advances  in  the  technologies  of  manufacturing  and 
materials  have  enhanced  the  current  application  of  composite  materials  from  being  used 
as  secondary  structural  elements  to  becoming  primary  load-carrying  structural 
components.  Due  to  the  inherent  inhomogeneity  and  anisotropy  of  these  materials, 
analysis  of  composite  structures  imposes  new  challenges  on  engineers. 

Plate  and  shell  structures  made  of  laminated  composite  materials  have  often  been 
modeled  as  an  equivalent  single  layer  using  classical  laminate  theory  (C.L.T.) 
[Christensen,  1979;  Jones,  1975;  Pagano,  1989  &  1973;  Pipes  &  Pagano,  1970]  in  which 
thickness  stress  components  are  ignored.  C.L.T.  is  a  direct  extension  of  classical  plate 
theory  in  which  the  well-known  Kirchhoff-Love  kinematic  hypothesis  is  assumed 
enforced.  This  theory  is  adequate  when  the  thickness  (to  side  or  radius  ratio)  is  small. 
However,  laminated  plates  and  shells  made  of  advanced  filamentary  composite  materials 
are  susceptible  to  thickness  effects,  because  their  effective  transverse  moduli  are 
significantly  small  compared  to  the  effective  elastic  modulus  along  the  fiber  direction. 
Furthermore,  the  classical  theory  of  plates,  which  assumes  that  the  normals  to  the  mid¬ 
plane  before  deformation  remain  straight  and  normal  to  the  plane  after  deformation, 
under  predicts  deflections  and  over  predicts  natural  frequencies  and  buckling  loads. 
These  discrepancies  arise  due  to  the  neglect  of  transverse  shear  strains.  The  range  of 
applicability  of  the  C.L.T.  solution  has  been  well  established  for  laminated  flat  plats 
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[Noor  &  Burton,  1989;  Pagano,  1989;  Reddy,  1993].  In  order  to  overcome  the 
deficiencies  in  C.L.T.,  refined  laminate  theories  have  been  proposed  [Byun  &  Kapania, 
1992;  Noor  &  Burton,  1990  &  1989;  Reddy,  1993,  Oct  and  Dec  1984;  Robbins,  Reddy  & 
Murty,  1991;  Whitney  &  Pagano,  1970].  These  are  single  layer  theories  in  which  the 
transverse  shear  stresses  are  taken  into  account.  They  provide  improved  global  response 
estimates  for  deflections,  vibration  fi'equencies  and  buckling  loads  of  moderately  thick 
composites  when  compared  to  the  classical  laminate  theory.  A  Mindlin  type  first-order 
transverse  shear  deformation  theory  (S.D.T.)  was  first  developed  by  AVhitney  and  Pagano 
[1970]  for  multi-layered  anisotropic  plates,  and  by  Dong  and  Tso  [1972]  for  multi¬ 
layered  anisotropic  shells.  Both  approaches  (C.L.T.  and  S.D.T.)  consider  all  layers  as 
one  equivalent  single  anisotropic  layer,  thus  they  can  not  model  the  warping  of  cross- 
sections,  that  is,  the  in-plane  distortion  of  the  deformed  normal  due  to  transverse  shear 
stresses.  Furthermore,  the  assumption  of  a  non-deformable  normal  results  in 
incompatible  shearing  stresses  between  adjacent  layers.  The  later  approach  also  requires 
the  introduction  of  an  arbitrary  shear  correction  factor,  which  depends  on  the  lamination 
parameters  for  obtaining  accurate  results.  It  is  well  established  that  such  a  theory  is 
adequate  to  predict  only  the  gross  behavior  of  laminates. 

The  exact  analyses  performed  by  Pagano  [1989;  1973;  Pipes  &  Pagano,  1970]  on  the 
composite  flat  plates  have  indicated  that  the  in-plate  distortion  of  the  deformed  normal 
depends  not  only  on  the  laminate  thickness,  but  also  on  the  orientation  and  the  degree  of 
orthotropy  of  the  individual  layers.  Therefore  the  hypothesis  of  non-deformable  normals, 
while  acceptable  for  isotropic  plates  and  shells,  is  often  quite  unacceptable  for  multi¬ 
layered  anisotropic  plates  and  shells  that  have  a  large  ratio  of  Young’s  modulus  to  shear 
modulus,  even  if  they  are  relatively  thin  [F.  Chang,  Perez,  &  K.Y.  Chang,  1990;  Reddy, 
1993,  Oct  &  Dec  1984].  Thus  a  transverse  shear  deformation  theory  which  also 
accounts  for  the  warping  of  the  deformed  normal  is  required  for  accurate  prediction  of 
the  elastic  behavior  (deflections,  thickness  distribution  of  the  in-plane  displacements, 
natural  fi’equencies,  etc.)  of  multi-layered  anisotropic  plates  and  shells. 

In  view  of  these  issues  we  have  developed  a  multi-director  shell  theory  for  composite 
laminates  with  the  following  attributes: 

1.  The  displacement  field  is  continuous  through  the  thickness  of  the  multi-layer 
structure  while  the  rotation  field  is  layer-wise  continuous  (in  2-D)  and  can  be 
discontinuous  across  the  finite  element  layers  through  the  thickness  direction.  This 
displacement  field  fulfills  a  priori  the  geometric  continuity  conditions  between 
contiguous  layers. 

2.  The  assumed  displacement  field  is  capable  of  modeling  the  distortion  of  the  deformed 
normal,  without  increasing  the  order  of  the  partial  differential  equations  with  respect 
to  the  first-order  transverse  shear  deformation  theory. 

3.  The  assumed  displacement  field  has  a  3-D  feature,  thereby  modeling  accmately  the 
interlaminar  conditions  and  predicting  the  3-D  edge  effects,  in  orthotropic  layers. 
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4.  Like  the  shear  deformable  theory,  the  proposed  composite  shell  theory  provides 
flexibility  in  the  specification  of  the  boundary  conditions. 

5.  In  this  theory,  at  most,  only  first  derivatives  of  displacement  and  rotation  fields 
appear  in  the  variational  equations.  The  practical  consequence  of  this  fact  is  that  only 
C®  continuity  of  finite  element  functions  is  required,  which  is  readily  satisfied  by  the 
family  of  LaGrange  elements. 

6.  Development  of  electromechanical  and  thermomechanical  models  for  simulation  of 
active  surfaces  with  sensors  and  actuators  elements  (either  embedded  or  mounted)  is 
straight  forward  and  it  is  a  natural  extension  of  the  theory  for  coupled  systems. 

Considerable  attention  has  been  devoted,  during  the  past  few  years  to  the  developments 
of  active  structural  systems  with  shape  control  capabilities  [Pfaeffle,  1993],  [Beauchamp, 
1992],  [McLean,  1993].  These  structures  range  fi-om  simple  beams  that  are  controlled 
with  a  single  actuator  to  the  more  imaginative  compliant  fins  that  are  controlled  with  a 
distributed  network  of  actuators.  In  these  structures,  the  emphasis  has  been  placed  on 
utilizing  Nitinol  fibers  which  are  embedded  either  directly  or  indirectly  inside  the  fabric 
of  these  multilayer  structures.  Generally  Nitinol  fibers  are  thermally  trained  to  shrink  and 
remain  straight  upon  heating  above  their  austenite  phase  transformation  temperature. 
Restraining  the  fibers  from  shrinking,  by  the  composite  when  the  fibers  are  directly 
embedded  or  by  end  restraints  when  embedded  indirectly  results  in  the  generation  of 
large  forces.  By  virtue  of  the  spatial  spacing  between  the  direction  of  the  developed 
phase  recovery  forces  and  the  neutral  planes,  control  moments  are  generated  which  are 
then  used  to  control  the  shape  of  the  structures.  In  this  manner,  the  shape  memory  effect 
is  not  fully  utilized  to  its  complete  potential  particularly  because  the  deflection  of  the 
Nitinol  fibers  is  limited  to  motions  along  the  fibers’  longitudinal  axes.  Furthermore,  the 
use  of  directly  embedded  fibers  results  in  thermal  buckling  of  the  composites  due  to  the 
thermal  stresses  induced  by  the  activation  of  the  Nitinol  fibers.  In  the  case  of  the 
indirectly  embedded  fibers,  mechanical  as  well  as  thermal  buckling  are  developed  due  to 
the  generated  in-plane  recovery  forces.  In  either  case  and  before  the  occurrence  of 
buckling,  the  stiffiiess  as  well  as  the  natural  frequencies  of  the  Nitinol-reinforced 
composites  is  reduced  considerably. 

An  alternative  way  of  controlling  the  shape  of  Nitinol-reinforced  composites  relies  on  the 
operation  of  Nitinol  strips  which  are  indirectly  embedded  inside  the  composites  [Baz, 
1994],  where  a  full  utilization  of  the  shape  memory  effect  to  control  the  shape  of  beams 
without  compromising  their  structural  stiffiiesses  or  fi-equencies  is  demonstrated. 

Several  biomedical  applications  have  also  resulted  in  innovative  designs  where  active 
structure  has  experienced  deformations  that  are  an  order  of  magnitude  larger  than  the 
actuator  material  [Barrett,  1995]. 
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To  develop  modeling  capabilities  suitable  for  these  applications,  a  comprehensive 
literature  review  was  conducted  on  geometrically  nonlinear  finite  element  approaches  and 
nonlinear  solution  techniques  proposed  in  the  last  two  decades.  Shape  control 
applications  are  mostly  associated  with  geometrically  nonlinear  kinematics  but  small 
strain  values.  As  a  result,  an  element  co-rotational  procedure  was  selected  along  with  a 
nonlinear  solution  technique  for  handling  the  geometric  nonlinearities. 

The  theoretical  issues  related  to  the  implementation  of  the  co-rotational  procedure  and  the 
close  relationship  between  the  Updated  Lagrangian  formulation  and  the  co-rotational 
procedure  are  discussed  in  Part  III  of  this  report.  All  the  approximations  made  in  the  co- 
rotational  procedure  as  well  as  the  restrictions  and  limitations  caused  by  them  are  also 
identified. 


4.0  Summary 

Part  II  of  this  report  presents  the  constitutive  modeling  of  shape  memory  alloys  and  its 
numerical  implementation.  Martensitic  phase  transformations  are  discussed.  Evolution 
equations  for  production  of  austensite,  and  single  variant/multiple  variant  martensite  are 
presented.  State  determination  for  finite  element  modeling  of  the  theory  is  discussed. 
Truss  and  beam  elements  of  SMA  for  both  linear  and  nonlinear  kinematics  are  developed 
and  a  series  of  numerical  simulations  are  performed. 

Part  III  presents  the  finite  element  formulation  of  geometrically  nonlinear  multi-layered 
shells  in  a  co-rotational  fi-amework.  Total  and  Updated  Lagrangian  formulations  for  the 
element  co-rotation  procedure  are  discussed.  Stabilization  techniques  for  under¬ 
integrated  nonlinear  shell  elements  are  investigated,  four-node  quadrilateral  and  eight- 
node  hexahedral  elements  are  developed.  Several  simulations  of  flat  and  curved 
structures  imder  various  loadings  are  performed  and  the  results  are  compared  with 
benchmark  problems  to  show  the  accuracy  of  the  proposed  theory. 

Publications  and  presentations  associated  with  this  work  are  listed  below. 

1.  “Finite  Strain  Formulation  of  Pseudoelastic  Materials,”  Computer  Methods  in 
Applied  Mechanics  and  Engineering,  148  (1997)  23-37. 

2.  “Multi-layered  Shell  Formulation  of  Composite  Laminates,”  International  Journal  for 
Numerical  Methods  and  Engineering  (in  press). 

3.  “Coupled  Thermomechanical  Simulation  of  Shape  Memory  Alloys,”  Smart  Structures 
and  Materials  Conference,  San  Diego,  March  1997. 

4.  “Finite  Element  Implementation  of  Large  Deformation  Coupled  Thermomechanical 
Response  of  Shape  Memory  Alloys.”  Fourth  US  National  Congress  on 
Computational  Mechanics,  San  Francisco,  August  1997. 
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5.  “Thermomechanical  Simulation  of  Expansion  and  Contraction  of  Shape  Memory 
Alloy  Stent,”  1®'  International  conference  on  Advanced  Biomaterials,  October  1997, 
Canada. 
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Part  II 

Shape  Memory  Alloys 
Theory  and  Numerical 
Implementation 


Chapter  1 
Introduction 


Part  II  of  the  report  covers  our  research  in  the  area  of  constitutive  modeling  of  shape  memory 
alloys,  applications,  general  characteristics  of  martensitic  transformations  and  a  new  con¬ 
stitutive  model  for  solid-solid  phase  transformations.  A  new  algorithm  is  developed  for  the 
thermo-mechanical  constitutive  model  based  on  trial  values  of  stress  and  strain.  Lastly,  the 
model  is  implemented  into  a  multi-dimensional  truss  and  a  two-dimensional  beam  element. 
For  each  of  these  cases  we  consider  both  linear  and  finite  kinematics. 

The  objective  of  this  part  of  the  report  is  the  formulation  and  implementation  of  a  robust 
and  efficient  class  of  finite  elements  for  the  treatment  of  first  order  solid  state  martensitic 
transformations.  A  constitutive  model  describing  the  macroscopic  behavior  of  martensitic 
transformations  is  developed  in  the  setting  of  linear  and  finite  kinematics.  Due  to  its  broad 
use  in  industry,  Nickel-Titanium  (NiTi)  is  considered  for  the  majority  of  the  numerical 
simulations.  However,  it  should  be  noted  that  these  developments  are  applicable  to  other 
binary  and  ternary  alloys. 

The  remainder  of  this  chapter  describes  the  motivation  of  the  present  work  as  well  as  a  brief 
literature  review  of  martensitic  transformations  from  a  material  and  numerical  standpoint. 
Lastly,  an  overview  of  the  subsequent  chapters  is  outlined. 
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1.1  Motivation 


For  the  past  several  decades,  researchers  have  been  working  in  the  area  of  solid-solid  phase 
transformations,  specifically  martensitic  transformations.  Some  alloys  which  undergo  marten¬ 
sitic  phase  transformations  have  the  ability  to  deform  up  to  a  strain  of  approximately  11% 
for  single  crystals  and  approximately  8%  for  polycrystals  Funakubo  [1984]  and  Duerig  et.al. 
[1990],  which  is  an  order  of  magnitude  greater  than  common  engineering  alloys  used  for 
design  structures  such  as  carbon  and  stainless  steels.  These  alloys  are  commonly  referred 
to  as  Shape  Memory  Alloys  (SMA).  The  range  of  strain  for  SMA,  although  moderate,  can 
be  “fully”  recovered  under  either  a  load  and/or  a  thermal  cycle.  The  conditions  for  this 
recovery  depend  upon  a  variety  of  parameters,  to  be  discuss  in  detail  in  Chapter  2. 

During  the  transformation  process,  the  crystal  structure  of  the  material  is  altered,  result¬ 
ing  in  microscopically  and  macroscopically  measurable  changes  in  the  mechanical,  thermal, 
acoustical,  electrical  and  magnetic  properties  Duerig  et.al.  [1990].  These  characteristics, 
which  can  be  designed  into  the  material  by  adding  alloying  elements,  make  SMA  attractive 
from  an  engineering  standpoint.  In  addition  to  the  physical  changes,  the  alloys  have  other 
features  which  make  them  attractive,  such  as  biocompatibility,  high  resistance  to  corrosion 
and  fatigue,  and  a  high  force  to  weight  ratio. 

There  are  a  variety  of  SMA  currently  available,  but  few  exhibit  the  range  of  motion  described 
above,  especially  under  cyclic  conditions.  The  most  common  alloys  used  are  a  binary  alloy 
Nickel-Titanium  (NiTi)  and  ternary  copper  alloys  CuZnAl  and  CuAlNi.  The  alloy  which  is  of 
most  interest  for  the  current  study  is  NiTi  commonly  referred  to  as  Nitinol  (Nickel  Titanium 
Naval  Ordnance  Laboratory),  due  to  its  high  recoverable  strain,  manufacturability  and  its 
potential  in  the  medical  industry. 

As  an  example  of  a  simple  engineering  application  for  SMA,  we  consider  the  area  of  smart 
structures,  namely  active  vibration  control  and  shape  control.  Through  a  variation  in  the 
kinetic  or  thermal  fields,  the  strain  level  and  material  properties  can  be  altered  to  facilitate 
a  desired  response  of  the  structure.  Figure  1.1  depicts  the  use  of  these  alloys  or  “smart 
materials”  in  a  simple  mass-spring  system  for  vibration  control. 

Figure  1.1  depicts  a  spring-mass  system  in  two  different  temperature  states.  The  system  on 
the  left  is  at  a  low  temperature,  T  =  Tq  where  the  crystal  structure  of  the  spring  is  such 
that  the  stiffness  is  given  as  AT  =  Kq.  The  system  on  the  right  is  at  a  high  temperature, 
T  =  To  +  AT,  under  the  same  state  of  stress.  At  the  elevated  temperature  the  crystal 
structure  of  the  spring  is  alter  such  that  the  resulting  stiffness  is  A"  =  aKo.  The  fundamental 
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Figure  1.1:  Vibration  control  for  a  spring-mass  system  using  smart  materials  for  a  constant 
load  and  variable  temperature. 

frequency  of  this  simple  spring-mass  system  is  given  by  w  =  \JK/M.  We  see  that  ratio  in 
the  fundamental  frequency  between  the  two  states  is  y/a.  As  an  example,  relative  changes 
in  the  frequency  as  high  as  a  factor  of  2  (i.e.  o;  w  4)  can  occur  for  temperature  variations  of 
AT  =  15  —  40(7  for  NiTi,  Funakubo  [1984]. 

As  previously  mentioned,  the  potential  for  SMA  in  the  medical  industry  is  increasing.  Since 
biocompatibility  is  paramount  for  medical  devices,  the  composition  of  the  alloys  must  be 
chosen  judiciously.  Fortunately,  one  of  the  best  biocompatible  alloys  in  the  class  we  are 
concerned  with  is  NiTi.  The  main  concern  with  regard  to  biocompatibility  or  toxicity  is 
the  surface  of  the  instrument  or  device  and  since  NiTi  naturally  forms  an  oxidized  layer  of 
Ti02  on  its  exterior,  it  is  a  good  candidate  for  medical  applications.  One  such  application 
is  a  medical  device  referred  to  as  a  stent,  see  Figure  1.2.  The  stent  is  placed  within  an 
artery  typically  after  full  balloon  angioplasty  to  aid  and  give  support  to  deteriorated  regions. 
Traditionally,  in  this  procedure  a  stainless  steel  stent  structure  is  placed  over  the  balloon 
and  expanded  with  the  balloon  against  the  artery  wall,  once  the  balloon  is  deflated  and 
removed  the  stent  remains  expanded.  Although  similar  in  design  and  construction  to  their 
stainless  steel  counterparts,  the  NiTi  stents  give  the  advantages  of:  deliverability  to  the 
region  of  choice  after  the  angioplasty  has  been  performed,  a  reduction  in  the  pressure  the 
stent  exerts  on  the  vessel  wall  (typically  the  stent  is  design  such  that  the  resulting  stiffness  of 
the  in-situ  stent  is  low),  it  detectability  under  exploratory  procedures  such  as  MRI,  increased 
biocompatibility,  the  ability  to  recover  its  geometry  in  the  event  of  a  local  bifurcation  of  the 
artery  due  to  some  external  stimuli  and  its  amenability  to  surface  coatings  for  local  drug 
delivery. 

Further  applications  are  discussed  in  Wayman  [1980]  and  Funakubo  [1984],  which  include 
package  clamps,  disc  seals,  shut  off  valves,  orthodontic  arch  wires,  medical  guidewires  and 
catheters,  eyeglass  frames,  orthopedic  implants  and  more. 

The  key  focus  of  the  present  work  is  to  develop  a  robust  constitutive  model,  which  is  readily 
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Figure  1.2:  NiTi  medical  stent. 


implemented  into  a  numerical  setting  for  the  analysis  and  design  of  shape  memory  alloy 
structures.  The  particular  avenue  used  to  numerically  solve  the  representative  equations  of 
the  model  problem  will  be  the  finite  element  method. 


1.2  Background 


The  first  recorded  observation  of  a  shape  memory  alloy  was  in  1932  by  Chang  &  Read  [1951] 
for  AuCd.  Subsequent  discoveries  of  various  copper  and  iron  based  alloys  were  discovered 
shortly  afterward.  However  in  1961  the  Naval  Ordnance  Laboratory  observed  that  NiTi  had 
properties  of  a  SMA,  Buehler  et.  al.  [1963].  Since  that  time  significant  work  has  been 
directed  to  the  understanding  of  both  the  microscopic  and  macroscopic  behavior  of  these 
alloys. 

There  are  a  variety  of  approaches  researchers  have  taken  to  model  and  understand  the 
behavior  of  shape  memory  alloys.  We  overview  some  of  the  more  outstanding  areas  and 
highlight  their  strategies. 

One  approach  is  concerned  with  describing  the  process  of  the  phase  transformation  on  a 
local  level.  This  approach  has  been  to  develop  models  based  upon  non-convex  multi-well 
free  energy  formulations  with  an  addition  of  a  non-local  interaction  terms  to  account  for 
interfacial  energy  at  the  phase  boundary.  Using  a  Landau-Devonshire  based  model  Falk 
[1980]  was  the  first  to  apply  a  multi-well  free  energy  function  to  describe  both  stress  and 
thermally  induced  phase  transformations  for  martensitic  transformations.  Additional  work 
soon  followed  along  these  lines  and  are  discussed  in  Falk  [1984],  Ball  &  James  [1992],  Sun 
&  Hwang  [1993ab],  Abeyaratne  &  Knowles  [1993],  Kafka  [1994ab],  Abeyaratne  et.al.  [1994], 
Kim  &  Abeyaratne  [1995],  Rogers  [1996],  and  Patoor  et.  al.  [1996].  Muller  &  Xu  [1991] 
extended  the  developments  above  to  model  the  characteristics  of  the  phase  transformations 


at  high  temperatures,  specifically  pseudoelasticity. 

Utilizing  arguments  from  thermodynamics  and  statistical  mechanics  Muller  &  Wilmanski 
[1981],  Achenbach  &  Muller  [1985]  Achenbach  [1986]  and  Achenbach  [1989]  have  developed 
a  rate-dependent  constitutive  model  to  describe  the  martensitic  phase  transformation  under 
kinetic  and  thermal  loading. 

The  significant  work  done  in  the  area  of  micro-mechanics  and  statistical  thermodynamics 
with  respect  to  martensitic  transformations  has  enabled  the  development  of  macro-structural 
models.  Specifically,  a  phenomenological  approach  in  which  the  micro-structure  is  accounted 
for  by  the  introduction  of  internal  variables  into  the  constitutive  theory.  Tanaka  k  Iwasakai 
[1985]  AND  Tanaka  [1986]  incorporated  the  earlier  work  on  free  energy  functions  to  develop  a 
rate-independent  model  which  describes  the  phenomenological  behavior  of  martensitic  phase 
transformations  for  a  one-dimensional  bar  under  a  tensile  stress  state.  Liang  k  Rogers  [1990] 
extending  this  work  by  integrating  Tanaka’s  rate  equations  to  develop  a  model  in  which 
evolution  equations  have  an  explicit  form.  Brinson  k  Lammering  [1993]  incorporated  the 
Tanaka  model  into  a  finite  element  setting  with  finite  kinematics  for  a  one-dimensional 
bar,  again  restricted  to  a  tensile  state  of  stress.  Ivshin  k  Pence  [1994ab]  extended  the 
constitutive  theory  to  incorporate  thermal  effects  under  a  tensile  state  of  stress,  but  due  to 
the  complexity  and  parameters  needed  for  their  model  (i.e.  specific  entropies),  their  model 
has  not  attracted  attention  from  an  implementation  standpoint.  A  general  phenomenological 
model  for  SMA  was  also  developed  by  Auricchio  [1995]  for  both  tensile  and  compressive  states 
for  pseudoelasticity. 

As  seen  above,  efforts  have  just  recently  begun  in  the  area  of  computational  mechanics  in 
regard  to  general  phenomenological  shape  memory  alloys  models.  The  numerical  models 
which  have  been  developed  thus  far  are  limited  to  a  high  temperature  regimes  and  account 
for  micro-structure  production/depletion  in  tensile  states  of  stress  only.  To  overcome  these 
deficiencies  we  have  developed  a  one-dimensional  constitutive  model  and  algorithm  which 
encompasses  the  fully  thermomechanical  regime  of  the  phase  transformation  space  for  both 
tensile  and  compressive  states  of  stress.  The  constitutive  model  is  based  upon  the  intro¬ 
duction  of  internal  variables  for  both  tensile  and  compressive  states  from  which  we  develop 
evolution  equations  which  capture  the  phenomenological  behavior  of  the  martensite  phase 
transformations  for  all  temperature  and  stress  states.  In  addition,  an  algorithm  for  state 
determination  is  developed  which  relies  on  modified  trial  state  variables. 
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1.3  Applications 


Shape  memory  alloy  devices  have  been  used  commercially  for  over  twenty  years,  Wayman 
[1980].  With  the  advent  of  shape  memory  alloys  in  the  market  place,  products  have  exhibited 
enhanced  functionality  with  respect  to  their  ability  to  react  and  adjust  to  an  external  process 
in  real  time.  These  products  have  a  wide  range  of  applications.  The  two  main  areas  of  focus 
are:  commercial  and  medical  applications. 

There  are  two  main  mechanisms  employed  for  the  design  and  construction  of  SMA  devices, 
free  and  constrained  recovery.  In  some  designs  a  combination  of  both  is  used.  Free  recovery 
typically  consists  of  an  initial  configuration  in  which  the  specimen  is  crystallographically 
martensite  and  under  the  action  of  external  stimuli  is  deformed  into  an  intermediate  con¬ 
figuration.  The  specimen  is  then  heated  and  freely  recovers  the  imposed  deformation,  at 
the  new  configuration  the  specimen  is  in  an  austenitic  state.  This  process  may  be  repeated 
by  cooling  the  specimen  from  the  austenitic  state,  at  zero  load,  to  form  martensite,  which 
is  then  loaded  and  subsequently  heated.  This  mode  of  recovery  is  widely  used  in  medical 
instruments  for  the  recovery  of  blood  clots  or  other  inclusions  within  the  body.  Constrained 
recovery  consists  again  of  an  initial  configuration  in  which  the  specimen  is  crystallographi¬ 
cally  martensite  and  under  the  action  of  external  stimuli  is  deformed  into  an  intermediate 
configuration.  The  specimen  is  then  placed  in  contact  with  another  component  in  the  sys¬ 
tem.  Upon  heating,  the  specimen  attempts  to  recover  the  deformation  initially  imposed  on 
it,  but  the  contact  points  restrict  a  full  recovery.  The  resulting  contact  points  firmly  join 
the  two  parts  together.  This  approach  is  widely  used  for  couplers  and  clamping  devices. 


1.3.1  Commercial  Devices 

Commercial  applications  include  some  of  the  following  products:  brassiere  underwires,  an¬ 
tenna  rods  for  cellular  telephones,  orthodontic  bridge  wires,  temple  and  bridge  components 
for  eye  glass  frames,  pipe  couplings,  coffee  maker  components,  thermostatic  mixing  valves, 
active  vibration  and  shape  control  of  composite  structures  embedded  with  shape  memory 
alloy  fibers  and  deployment  and  release  mechanisms  in  extreme  conditions,  such  as  space 
and  deep  sea  exploration. 

The  first  significant  application  in  the  commercial  sector  were  pipe  couplers,  Wayman  [1980] 
employing  constrained  recovery  to  join  two  pipes  or  repair  cracks  in  existing  pipelines, 
see  Figure  1.3  The  process  for  the  coupling  begins  initially  at  ambient  temperature  in  an 
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austenitic  state.  The  coupling  is  cooled  below  a  critical  value  and  loaded  and/or  deformed 
mechanically  until  the  coupling  has  transformed  into  100%  martensite  crystallographically. 
The  expanded  coupling  is  placed  over  the  location  to  be  joined  and  allowed  to  heat  up  to 
ambient  temperature.  During  the  thermal  cycle  a  reverse  phase  transformation  to  austenite 
occurs,  along  with  the  associated  strain  recovery.  The  resulting  parts  are  securely  fasten  by 
the  SMA  coupling.  Years  of  service  have  shown  remarkable  reliability  of  these  shape  memory 
alloy  connectors  under  various  conditions  leading  to  the  extended  use  of  SMA  constrained 
recovery  for  other  sealing  and  fastening  devices. 


Figure  1.3;  Shape  memory  alloy  pipe  coupling.  The  process  begins  with  the  specimen  at 
ambient  temperature  (1)  and  is  then  cooled  at  a  constant  load  (typically  zero)  transforming 
the  austenite  to  multiple  variant  martensite  (2).  The  coupling  is  then  expanded,  transforming 
the  specimen  to  a  single  variant  of  martensite  and  placed  over  the  pipe  to  be  joined  (3). 
Lastly,  the  coupling  is  allow  to  heat  to  ambient  temperature  and  contracts  (4). 


1.3.2  Medical  Devices 

The  medical  applications  are  numerous  and  has  potential  for  significant  growth  in  appli¬ 
cations.  The  predominant  shape  memory  alloy  structures  used  in  the  medical  field  are 
guidewires,  hingeless  and  steerable  instruments  and  self  expanding  stents. 

Guidewires  are  used  to  access  the  desired  location  in  the  body  from  an  external  entry  point. 
Once  the  guidewire  is  in  place,  the  device,  such  as  a  catheter,  is  delivered  over  the  guidewire  to 
area  of  interest.  With  its  inception  in  1953  (Fernald  et.  al.  [1994]),  stainless  steel  guidewire 
have  aided  the  surgeon  in  many  tasks.  The  recent  introduction  of  Nitinol  guidewires  in 
1987  (Fernald  et  al  [1994])  has  gained  wide  acceptance  due  to  its  superior  performance 
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over  its  stainless  steel  counterpart.  The  difference  of  shape  memory  alloy  guidewire  over 
typically  used  stainless  steel  is  the  superior  performance  in  trackability,  kink  resistance, 
flexibility  and  torquability.  These  features  yield  a  guidewire  which  can  maneuver  through 
various  pathways  in  arteries  without  undergoing  permanent  deformation  or  buckling,  while 
maintaining  excellent  controllability. 

Self-expanding  stents,  see  Figure  1.2,  are  utilized  for  the  treatment  of  stenosis  of  hollow 
organs  or  arteries  by  aiding  in  the  support  of  the  deteriorated  region,  see  Ozaki  [1996]  and 
SiGWART  [1997]  for  a  general  review  of  stents  and  new  technologies.  The  placement  of  the 
stent  is  described  below.  The  stent  with  an  initial  diameter  of  3  —  15mm  is  cooled  below  a 
critical  temperature  and  placed  into  a  delivery  tube  with  an  outside  diameter  of  1  —  6mm. 
The  delivery  tube  with  the  stent  inside  is  then  heated  to  room  temperature.  The  delivery 
tube  is  placed  on  a  guidewire  and  brought  to  its  desired  location.  Once  the  delivery  tube  is 
at  the  desired  location,  the  stent  is  pushed  out  of  the  delivery  system.  The  stent  undergoes  a 
reverse  transformation  (martensite  to  austenite)  and  expands  to  its  original  shape  and  size. 
The  stents  are  designed  such  that  the  final  configuration  minimizes  the  contact  pressure  with 
the  vessel  wall.  The  placement  of  the  stents  keeps  the  vessel  open  and  allows  unimpeded 
blood  flow.  The  advantage  over  the  stainless  steel  stent  model  is  its  ability  to  expand  and/or 
contract  twice  or  more  its  size  without  the  onset  of  permanent  deformation.  The  ability  of 
contracting  the  stent  without  permanent  damage  reduces  the  size  of  the  delivery  tube  in 
which  the  stent  is  placed,  which  in  turn  avoids  large  incisions  or  surgery  for  the  placement  of 
the  delivery  tube  and  reduced  trauma  to  tissue.  These  benefits  result  in  accelerated  recovery 
for  the  patient  and  reduced  hospital  costs. 


1.4  Overview 

Part  II  of  this  report  is  directed  toward  constitutive  theory  and  numerical  implementation 
using  the  finite  element  method  of  shape  memory  alloys.  An  outline  of  each  section  is  briefly 
discussed  below. 

In  Chapter  2  we  discuss  the  general  characteristics  of  martensitic  phase  transformations 
and  the  associated  macroscopic  behavior.  Chapter  3  summarizes  the  general  theory  of  a 
thermomechanical  continuum  with  the  inclusion  of  internal  variables.  Constitutive  models 
for  shape  memory  alloys  for  linear  and  finite  kinematics  are  also  developed  and  discussed 
in  Chapter  3.  Numerical  and  algorithmic  developments  with  respect  to  the  finite  element 
method  are  proposed  in  Chapter  4  for  a  multi-dimensional  truss-ba^  and  a  two-dimensional 
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beam  element.  Several  simulations  are  presented  which  show  the  qualitative  and  quantitative 
behavior  of  the  model  under  various  loading  conditions.  The  development  of  the  algorithm 
for  the  determination  of  state  is  also  presented,  along  with  its  advantages  over  commonly 
used  methods. 
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Chapter  2 

Martensitic  Phase  Transformations 


This  chapter  is  dedicated  to  describing  the  general  characteristics  of  martensitic  phase  trans¬ 
formations.  Both  mechanical  and  thermodynamic  aspects  of  phase  transformation  classifi¬ 
cation  and  activation  will  be  addressed.  We  also  discuss  the  characteristic  features  of  the 
shape  memory  effect  and  the  pseudoelastic  effect  for  shape  memory  alloys. 


2.1  General  Characteristics 


Martensitic  transformations,  classified  as  first  order  diffusionless  displacive  transformations, 
consist  of  lattice  transformations  involving  shearing  deformation,  which  results  from  a  co¬ 
operative  motion  of  atoms  over  a  small  distance.  This  movement  is  such  that  there  exists 
a  one-to-one  correspondence  between  lattice  points  in  the  parent  phase  (austenite)  to  the 
lattice  points  in  the  product  phase  (martensite),  known  as  lattice  correspondence.  The  inter¬ 
face  between  the  parent  and  the  product  phases  corresponds  to  the  plane  on  which  shearing 
occurs  during  the  transformation.  This  plane  is  commonly  referred  to  as  the  habit  plane  and 
during  transformation  is  accompanied  by  a  macroscopic  change  in  shape.  Throughout  the 
transformation,  the  habit  plane  experiences  no  strain  or  rotation.  These  conditions  lead  to 
preservation  of  planes  and  lines  between  the  parent  and  product  phases  and  is  thus  termed 
an  invariant  plane  which  can  be  described  by  a  linear  transformation  as  shown  in  Figure  2.1. 

During  a  martensite  transformation,  the  motion  shown  in  Figure  2.1  represents  an  ideal¬ 
ization.  Typically,  a  number  of  martensite  crystals  are  formed  with  differing  habit  plane 
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Figure  2.1:  Martensitic  transformation  from  a  multiple  variant  state  to  single  variant  state. 
The  solid  and  hatched  regions  represent  different  variants  of  martensite,  which  coalesce  into 
a  single  variant  under  external  stimuli.  The  plane  or  interface  between  the  variants  is  called 
the  habit  plane. 

indices,  but  all  crystallographically  equivalent.  The  differing  habit  planes’  indices  are  called 
variants. 

Also  note  from  Figure  2.1  the  martensitic  transformation  is  predominantly  a  shearing  type 
of  motion,  resulting  in  a  quasi-isochoric  motion.  The  relatively  small  change  in  volume 
(—0.1%  to  0.1%)  Fisher  [1994]  and  lattice  structure  allows  the  transformation  to  proceed 
without  incurring  plastic  deformation  in  the  parent  phase.  Whereas,  for  ferrous  alloys  the 
volume  change  is  typically  an  order  of  magnitude  larger  (although  regarded  as  being  iso- 
choric)  and  plastic  deformation  is  evident  in  the  parent  phase. 

Lastly,  the  crystal  structure  of  martensite  is  relatively  less  symmetric  compared  to  the  parent 
phase  from  which  various  variants  of  martensite  occur.  Due  to  the  reduced  symmetry,  the 
variety  of  lattice  correspondences  between  phases  involved  in  the  reverse  transformation 
are  restricted.  A  result  of  its  ordered  structure,  the  orientation  of  the  parent  phase  is 
automatically  preserved. 


2.2  Classification  of  Transformations 


There  are  two  types  of  martensitic  transformations:  athermal  and  isothermal.  For  the  devel¬ 
opment  of  the  shape  memory  alloy  model  we  only  need  to  consider  athermal  transformations. 
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outlined  below. 


An  athermal  transformation  is  a  transformation  which  advances  when  the  thermal  field 
is  below  some  critical  temperature  in  a  stress  free  reference  configuration.  Above 
this  temperature,  the  parent  phase  is  stable  (i.e.  the  material  has  a  natural  affinity  to  an 
austenitic  state).  The  transformation  may  be  non-thermoelastic  or  thermoelastic.  If  non¬ 
thermoelastic  during  the  transformation,  martensite  micro-regions  are  nucleated  and  grow 
to  their  final  size  over  a  very  small  time  interval  (approximately  1/3  the  elastic  wave  speed  of 
the  solid,  Fischer  [1994] ).  Continued  cooling  will  not  affect  further  growth  of  the  martensite 
crystals.  The  transformation  will  continue  by  further  nucleation  and  growth  of  martensite 
until  in  the  remaining  parent  phase  is  depleted.  Conversely,  in  thermoelastic  transformations 
after  the  martensite  micro-regions  are  nucleated  they  continue  to  grow  at  a  rate  proportional 
to  the  thermal  field.  Shape  memory  alloys  are  an  example  of  thermoelastic  transformations. 

Moreover,  most  thermoelastic  transformations  occur  in  ordered  alloys  since  they  are  crys- 
tallographically  reversible.  Since  the  original  orientation  of  the  parent  phase  is  necessarily 
selected  during  the  transformation,  alloys  which  undergo  thermoelastic  transformations  also 
experience  the  shape  memory  effect,  which  will  be  discussed  in  the  following  sections. 


2.3  Activation  of  Transformations 


Activation  of  a  martensitic  transformation  occurs  due  to  the  presence  of  a  driving  forces, 
either  thermal  of  kinetic.  To  initiate  a  transformation,  the  chemical  free  energy  difference 
between  the  parent  and  product  phases  must  be  greater  than  the  necessary  free  energy 
barriers,  such  as  transformational  strain  energy  or  interface  energy,  as  shown  in  Figure  2.2. 

The  temperature  T^s,  indicates  the  critical  value  the  temperature  must  be  reduced  to,  from 
the  equilibrium  temperature  Tq,  in  order  for  the  forward  transformation  to  occur.  This 
reduction  in  temperature  and  subsequent  value  for  the  free  energy  is  the  necessary  driving 
force  for  the  martensite  transformation.  The  temperature  denoted  by  Tas  indicates  the 
critical  temperature  at  which  the  reverse  transformation  will  occur,  if  the  temperature  is 
increased.  The  forward  and  reverse  transformation  can  be  seen  as  a  product  wedge  within 
the  parent  phase,  as  depicted  in  Figure  2.3,  under  the  action  of  an  external  driving  force 
either  thermal  or  kinetic. 

^The  symbol  T^g  denotes  the  critical  temperature  below  which  martensite  is  produced. 
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Figure  2.2:  Free  energy  versus  temperature.  At  To  we  see  that  both  austenite  and  martensite 
can  co-exist.  Whereas,  at  lower  temperatures  austenite  is  unstable  and  the  system  has 
a  natural  affinity  towards  martensite.  The  opposite  is  true  for  high  temperatures  where 
austenite  is  stable  and  martensite  is  not. 


Figure  2.3:  Product  wedge.  The  upper  diagram  shows  the  progression  of  the  martensite 
(arrow  shape  figures)  in  the  austenite  matrix  (solid  background)  under  the  action  of  external 
stimuli.  Whereas,  the  lower  diagram  shows  the  depletion  of  the  martensite  (arrow  shape 
figures)  in  the  austenite  matrix  (solid  background)  under  the  action  of  external  stimuli. 


For  either  the  forward  or  reverse  transformation  eventually  the  sum  of  the  chemical  and 
non-chemical  free  energies  approach  a  certain  minimum  value  and  growth  is  arrested  and 
the  transformation  is  complete. 

For  the  determination  of  when  transformations  initiate,  the  space  parameterized  by  stress 
and  temperature  is  commonly  used,  instead  of  the  space  parameterized  by  free  energy  and 
temperature.  The  stress-temperature  space  is  referred  to  as  the  phase  space  and  is  depicted 
in  Figure  2.4  for  tensile  states  of  stress  . 
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Figure  2.4:  Phase  space  diagram  depicting  transformation  zones  and  their  associated  values 
of  stress  and  temperature  for  activation. 


In  Figure  2.4  the  arrows  indicate  the  directionality  for  an  active  transformation,  while  cr®,. 
and  gIj.  denote  the  critical  start  and  finish  stress  of  a  martensite  transformation  and  Cm 
and  Ca  represent  of  the  slope  of  the  transformation  lines.  The  transformation  temperatures, 
Tmf,Tms,Tas  and  Taf  indicate  the  start  and  finish  temperature  at  zero  stress  for  martensite 
and  austenite,  respectively.  Lastly,  Gy  indicates  the  value  of  stress  above  which  plastic  slip 
will  occur.  Dependant  upon  the  path  taken  within  the  phase  space,  certain  characteristic 
features  of  the  stress-strain  space  (a  —  e)  will  be  manifested.  Further  discussion  of  various 
paths  is  addressed  in  the  following  sections. 


2.4  Shape  Memory  Effect 


If  a  martensite  transformation  is  induced  purely  by  a  thermal  field  below  a  critical  temper¬ 
ature,  Ta/^,  the  resulting  effect  is  called  “self  accommodation”  or  “twinning”.  The  ensuing 
multiple  variants  which  form  tend  to  average  the  overall  deformation  yielding  a  configuration 
similar  to  the  parent  configuration.  Macroscopic  deformation  is  not  observable  under  such 
a  transformation,  neglecting  the  thermal  expansion  term. 

If  the  martensite  transformation  is  further  induced  by  a  kinetic  field  the  multiple  variants 
which  are  present  will  coalesce  into  one  variant  in  the  preferred  direction  of  loading,  in  a 
process  known  as  detwinning.  Upon  removal  of  the  kinetic  field,  a  permanent  deforma- 

^The  symbol  To/  denotes  the  critical  temperature  above  which  austenite  is  stable. 
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tion  is  retained  until  the  specimen  is  heated  above  the  critical  temperature,  Taf,  leading 
to  a  recovery  of  the  residual  deformation.  Note,  to  obtain  a  full  reverse  transformation, 
the  martensite  transformation  must  undergo  a  thermoelastic  martensite  transformation  (i.e. 
crystallographically  reversible)  and  deformation  by  slip  must  hot  occur.  The  process  of  an 
initial  thermal/mechanical  loading  followed  by  a  thermal  loading  to  recover  the  transforma¬ 
tion  strain  is  typically  termed  the  shape  memory  effect.  The  basic  schematic  is  shown  in 
Figure  2.5. 


Figure  2.5:  Shape  memory  effect.  The  process  begins  stress  free  at  a  high  temperature  and 
is  cooled  under  zero  load  to  form  martensite.  The  specimen  is  then  loaded  and  unloaded  at 
constant  temperature.  Last,  the  specimen  is  heat  to  return  to  its  original  austenitic  state. 


2.5  Pseudoelastic  Effect 


When  the  thermal  field  is  above  the  critical  temperature  Taf  and  the  specimen  is  loaded 
mechanically  above  a  critical  stress  level  £7^/^,  the  austenite  crystal  will  transform  into  a 
single  variant  martensite  oriented  in  the  direction  of  loading,  accompanied  by  a  macroscopic 
strain  as  high  as  11%.  The  strain  is  recovered  upon  removal  of  the  mechanical  load,  since 
martensite  is  not  stable  at  low  stress  and  high  temperatures.  Typically  this  type  of  trans¬ 
formation  is  called  pseudoelasticity,  since  the  behavior  is  such  that  the  material  returns  to 

®The  symbol  amf  denotes  the  critical  stress  above  which  only  martensite  is  stable. 
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its  initial  configuration  upon  removal  of  the  loading.  The  basic  schematic  is  shown  in  Figure 

2.6 


Figure  2.6:  Pseudoelastic  process  begins  stress  free  at  a  high  temperature  and  is  then  loaded 
to  a  martensitic  state  and  unloaded  at  constant  temperature  to  its  original  austenitic  state. 

For  a  pseudoelastic  transformation  to  obtain  a  full  reverse  transformation  the  martensite 
transformation  must  undergo  a  thermoelastic  martensite  transformation  (i.e.  crystallograph- 
ically  reversible)  and  deformation  by  slip  must  not  occur. 


2-16 


Chapter  3 

One  Dimensional  Shape  Memory 
Alloy  Model 


In  this  chapter  we  develop  the  framework  for  a  solid-solid  phase  transformation  model. 
Internal  variables  are  used  to  phenomenologically  treat  the  response  of  the  microstructure 
on  a  macroscopic  level. 

We  begin  with  a  review  of  the  basic  equation  of  a  continuum.  Following  the  review,  linearized 
thermoelasticity  is  presented  in  one-dimension  for  subsequent  use.  The  phenomenological 
constitutive  equation  describing  the  solid-solid  phase  transformation  behavior  is  presented 
and  discussed.  To  complete  the  constitutive  equation  for  the  model  evolution  equations  for 
various  transformation  zones  within  the  phase  space  are  explicitly  derived. 


3.1  Continuum  Mechanics 


We  shall  give  a  brief  introduction  to  the  classical  theory  of  continuum  mechanics  in  the 
setting  of  nonlinear  elasticity.  For  a  more  in-depth  background  in  the  area,  an  excellent 
exposition  on  the  subject  of  continuum  mechanics  is  given  in  Truesdell  &  Noll  [1992],  as 
well  as  in  Gurtin  [1981].  This  section  is  given  for  completeness  and  is  used  as  a  framework 
for  later  use. 
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3.1.1  Kinematics 


Description  of  Motion 

Let  Ko^  denote  the  reference  configuration  of  a  deformable  bounded  body  B  CR^.  Identify 
a  material  point  X  with  its  place  X  in  the  configuration  by 


X  =  Ko{X)  . 


(3.1) 


Let  the  place  of  the  material  point  X  at  time  t  within  the  body  B  have  a  position  vector  x, 
given  by 


x  =  x{X,t),  (3.2) 

know  as  the  material  description  of  motion.  Using  (3.1)  and  (3.2)  and  assuming  invertibility 
of  Ko,  we  can  develop  a  relationship  between  the  place  of  material  point,  x  to  that  of  the 
material  point  at  time  t  as 


x  =  x«'(X),t)  =  x(X,t). 


(3.3) 


Thus,  the  material  point  X  is  mapped,  upon  deformation,  into  a  spatial  position  x  in  a 
current  configuration  at  the  current  time  by  means  of  a  single- valued,  continuously  differen¬ 
tiable  mapping  function  called  the  motion.  We  shall  also  assume  that  the  mapping  X  to  x 
is  one-to-one  and  onto,  hence  it  possesses  a  unique  inverse  for  all  time  t,  which  describes  the 
inverse  motion 


X  =  x^(x,t).  (3.4) 

^In  general,  the  subscript  zero  may  not  correspond  to  t  =  0  and  Kq  need  only  be  an  admissible  configu¬ 
ration,  not  necessarily  one  occupied  by  the  body. 
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A  necessary  and  sufficient  condition  for  the  invertibility  of  (3.3)  to  exist  is  for  the  Jacobian 
determinant  of  the  mapping  from  X  to  x  not  vanish  in  an  e  —  neighborhood  about  any  point 
within  the  body  B 


#  0  .  (3.5) 

We  shall  also  impose  the  physical  restriction  that  J  >  0,  ensuring  that  any  two  material 
points  do  not  occupy  the  same  position  in  space  at  the  same  time.  The  motion  described 
in  (3.3)  is  traditionally  referred  to  as  the  Referential  or  Lagrangian  description  of  motion, 
where  X  and  t  are  the  independent  variables,  while  (3.4)  is  called  the  Spatial  or  Eulerian 
description  of  motion,  with  x  and  t  as  the  independent  variables. 


Deformation  Gradients 

From  the  motion  described  in  (3.3)  we  define  the  deformation  gradient  F  relative  to  X  from 
some  fixed  reference  configuration  Ko  to  a  current  configuration,  Kt  as 


F  = 


dx 
dX  ■ 


(3.6) 


We  impose  the  following  restrictions;  that  oo  >  J  =  det{F)  >  0  to  ensure  that  no  region  of 
volume  about  an  e  —  neighborhood  of  a  material  point  becomes  zero  or  infinite  upon  defor¬ 
mation  and  the  derivatives  of  are  continuous.  The  deformation  gradient  characterizes  the 
transformation  of  a  line  element  dX  at  the  place  X  e  Kq  upon  deformation  to  a  line  element 
dx  at  the  place  x  6  kj.  Similarly,  from  the  inverse  motion  (3.4)  we  define  the  gradient  f 
relative  to  x  as 


f  =  ^  ,  (3.7) 

again  with  the  restrictions  that  oo  >  J~^  =  det{f)  >  0  and  the  derivatives  of  are 
continuous. 
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If  we  assume  for  convenience  that  the  reference  and  current  coordinate  frames  coincide,  then 
the  displacement  vector  u  from  the  two  coordinate  frames  can  be  expressed  as 


u  =  X  -  X  .  (3.8) 

Thus,  the  deformation  gradient  and  its  inverse  may  be  cast  in  terms  of  displacements  as 

F  =  I+^  =  I  +  H  and  f  =  I-^  =  I-h,  (3.9) 

where  H  and  h  are  called  the  Lagrangian  and  Eulerian  displacement  gradients,  respectfully. 
Defined  as 


H  =  F-  I  and  h  =  I-  f  =  . 


(3.10) 


Deformation  and  Strain  Measures 

To  introduce  a  deformation  measure  we  will  examine  the  magnitude  of  a  differential  line 
segment 


ds^  =  dyi-dx=  (FdX)  •  (FdX)  =  (FN)  •  (Fd'N)dS^  ,  (3.11) 


ds^ 

d^ 


=  =  N  •  F^FN  =  N  •  CN  , 


(3.12) 


where  A  is  the  stretch  of  the  differential  line  element.  Note  from  (3.12)  that  C  is  a  second 
order,  symmetric,  positive  definite  tensor,  which  is  a  common  measure  of  deformation  in  the 
reference  configuration  called  the  right  Cauchy- Green  deformation  iensor  defined  as 
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C  =  F^F  . 


(3.13) 


Similar  arguments  can  be  used  to  obtain  the  symmetric  positive  definite  second  order  tensor 
b,  which  is  a  common  measure  of  deformation  in  the  current  configuration  called  the  left 
Cauchy-Green  deformation  tensor  defined  as 


b  =  FF^  .  (3.14) 

When  a  motion  is  rigid,  the  line  segment  magnitudes  dS  and  ds  are  equal  and  thus  the  motion 
is  distance  preserving.  Since  the  line  segments  dS  and  ds  are  equal  for  a  rigid  motion,  the 
two  deformation  measures  are  also  equal  and  unity  (i.e.  C  =  b  =  I).  We  would  expect  a 
strain  measure  to  vanish  for  such  a  motion.  Hence,  we  shall  introduce  strain  measures  for 
the  reference  and  current  configurations.  Using  the  difference  of  the  magnitudes  of  the  line 
segments  squared  as  our  measure  of  strain,  we  obtain 


ds^  -  dS^  =  (C  -  I)  dX  dX  =  2  E  dX  dX  , 


(3.15) 


where  E  is  known  as  the  Green  or  Lagrangian  strain  tensor.  Defined  as 


E  =  i  (C  -  I)  .  (3.16) 

Similarly,  we  can  define  a  strain  measure  in  the  current  configuration  known  as  the  Almansi 
or  Eulerian  strain  tensor  by 


e  =  5(I-b-‘). 


(3.17) 
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Rotation  and  Stretch  Tensors 


Since  the  deformation  gradient  F  is  non-singular  by  (3.5),  the  polar  decomposition  theorem 
states  that  there  exist  positive  symmetric  tensors  U  and  V,  and  a  proper  orthogonal  tensor 
R  such  that 


F  =  RU  =  VR  .  (3.18) 

Moreover,  the  right  and  left  polar  decompositions  are  unique.  Where  the  tensors  U  and  V 
are  called  the  right  and  left  stretch  tensors  respectively,  while  R  is  called  the  rotation  tensor. 
Using  (3.13),  (3.14)  and  (3.18)  we  can  relate  the  right  Cauchy-Green  deformation  tensor  C 
to  the  right  stretch  U  and  the  left  Cauchy-Green  deformation  tensor  b  to  the  left  stretch  V 

by 


C  =  and  b  =  . 


(3.19) 


We  may  also  relate  U  and  V  by 


V  =  RUR^  .  (3.20) 

Note,  as  a  consequence  of  the  similarity  transformation  of  (3.20),  U  and  V  have  the  same 
eigenvalues. 


3.1.2  Balance  Laws 

Conservation  of  Mass 

We  assume  there  exists  a  function  m(f2),  termed  the  mass  of  the  body  Q  such  that 


m(fi)  >  0 


(3.21) 
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where  m{Q,)  is  independent  of  the  configuration  Kt  of  the  body  such  that 


=  0  .  (3.22) 

at 

The  function  m(f2)  may  be  expressed  in  terms  of  a  density  function  (a  scalar)  field  p  defined 
for  a  particular  configuration  as 


m(Q)  =  /  p{x,t)dv  =  0  (3.23) 

Jnt 

where  x  is  the  place  occupied  by  the  material  point  X  in  the  configuration  Kt,  t  is  time, 
dv  is  the  volume  element  corresponding  to  Kt  and  p  is  the  mass  density  of  the  body  in  the 
configuration  kj,  assumed  to  be  smooth  and  continuous.  Since  m(f])  is  independent  of  the 
configuration  ,  we  may  express  equation  (3.23)  as 


dV 


(3.24) 


for  an  arbitrary  reference  configuration  kq)  where  X  is  the  place  occupied  by  the  material 
point  X  in  the  configuration  kq,  dV  is  the  volume  element  corresponding  to  kq  and  po  is  the 
mass  density  of  the  body  in  the  configuration  kq.  Noting  that  (3.24)  holds  for  an  arbitrary 
body  n  and  the  mass  density  is  smooth  and  continuous,  we  may  express  the  local  form  of 
(3.24)  as 


Jp  =  Po  (3.25) 

where  we  have  used  the  relation  dv  =  JdV  which  relates  a  material  volume  in  the  reference 
configuration  kq  to  an  arbitrary  configuration  «<,  J  given  in  the  previous  section.  We  may 
take  the  time  rate  of  change  of  (3.24)  to  yield 
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d  C 

— J  p{x,  t)  dv  =  0 


(3.26) 


or  locally  as 


p  +  /xiiv  [v]  =  0 


(3.27) 


where  the  superposed  dot  indicates  the  material  time  derivative.  Equations  (3.24)  through 
(3.27)  are  all  different  forms  of  the  conservation  of  mass. 


Linear  and  Angular  Momentum 

The  balance  of  linear  momentum  is  expressed  as  an  equality  between  the  rate  of  change  of 
linear  momentum  of  an  arbitrary  volume  of  the  body  to  the  resultant  force  acting  on  that 
volume  or  simply 


dt 


V  dm  = 


(3.28) 


where  dm  =  pdv  is  the  element  of  mass,  da  is  the  area  element,  v  is  the  velocity,  b  is  the 
body  force  per  unit  mass,  t  is  the  contact  force  per  unit  area,  V  is  an  arbitrary  volume  of 
the  body  B  and  dV  is  the  boundary  surface  of  V. 

Likewise,  the  balance  of  angular  momentum  is  expressed  as  the  balance  between  the  rate  of 
change  of  angular  momentum  for  an  arbitrary  volume  and  the  resultant  moment  of  the  force 
acting  on  that  volume  or 


dt 


X  X  V  dm  = 


X  X  b  dm  + 


X  X  t  da 


where  x  is  the  position  vector  from  the  origin  o. 


(3.29) 
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Recall  Cauchy’s  theorem:  There  exists  a  2"^^  order  tensor  field  T  independent  of  n  such  that 


t{x,n)  =  T{x)n  (3.30) 

for  all  X  E  13  and  for  arbitrary  unit  vectors  n.  The  tensor  T  is  referred  to  as  the  Cauchy 
stress  tensor. 

Substituting  (3.30)  into  (3.28)  and  utilizing  the  divergence  theorem,  the  continuity  of  p,  b 
and  V  and  assuming  T  is  once  continuously  differentiable  we  arrive  at 


div  [t^]  +  pb  =  pv  in  12  X  /  (3.31) 

where  Q.  is  the  domain  under  consideration,  div  [•]  is  the  divergence  operator  with  respect  to 
the  spatial  configuration  and  /  =  [0,  t]  represents  an  interval  of  time. 

Substitution  of  (3.27),  (3.30)  and  (3.31)  into  (3.29)  yields  symmetry  of  the  stress  tensor 


-T  in  X  /  .  (3.32) 

Upon  application  of  Nanson’s  formula^,  the  resultant  contact  forces  may  be  expressed  in 
terms  of  the  reference  configuration  and  subsequently  the  local  form  of  the  balance  of  linear 
momentum  is  expressed  as 


Div  [P]  +  pqB  =  poV  in  Qx  I  (3.33) 

where  P  is  the  first  Piola-Kirchhoff  stress,  B  is  the  body  force,  Div  [•]  is  the  divergence 
operator  with  respect  to  the  reference  configuration  and  V  =  k  is  the  material  velocity. 
While  the  local  form  of  the  balance  of  angular  momentum  is 


FP  =  p'^f'^  in  fix/. 

^Nanson’s  formula  relates  the  reference  and  current  area  elements  via  nda  =  JF~'^ NdA 


(3.34) 
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Balance  of  Energy 


The  balance  of  energy  or  the  first  law  of  thermodynamics  is  expressed  as  an  equality  between 
the  rate  of  change  of  the  internal  energy  plus  the  kinetic  energy  and  the  rate  of  working  of 
the  applied  forces  plus  any  heat  energy  entering  the  body 


A 

dt 


I  ( pe  +  ^pv  •  w  j  dv  =  j  {pr  +  pb  ■  v)  dv  +  j  {t  ■  v  —  h)  da  (3.35) 
JV  \  2  /  J-p  jQ-p 


where  p  is  the  current  mass  density,  e  is  the  specific  internal  energy,  v  is  the  velocity,  r 
specific  heat  supply,  b  is  the  body  force  per  unit  mass,  t  is  the  contact  force  per  unit  area, 
h  =  q  -  n\s  the  heat  flux  across  the  boundary  surface,  V  is  an  arbitrary  volume  of  the  body 
B  and  dV  is  the  boundary  surface  of  V. 

Utilizing  the  conservation  of  mass,  the  balance  of  linear  momentum,  Cauchy’s  theorem,  the 
transport  theorem  and  the  divergence  theorem  we  obtain 


/  {T  :  L  —  pe pr  —  div  [g])  dv  =  0  (3.36) 

Jr 

where  L  is  the  spatial  velocity  gradient. 

Since  (3.36)  holds  for  an  arbitrary  volume  V  of  the  body  B  and  we  assume  continuity  of  the 
arguments  within  the  integrand,  we  obtain  the  local  form  for  the  balance  of  energy 


T  :  L  —  pe-\-  pr  —  div  [g]  =  0  in  f2  x  /  (3.37) 

where  the  the  first  term  above  is  referred  to  as  the  stress  or  mechanical  power.  An  alternative 
material  form  for  the  balance  of  energy  expressed  in  terms  of  entropy  is  given  as 


{Tt])  =  -Div  [q]  +  R 


(3.38) 
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where  T  is  the  thermal  field,  tj  is  the  entropy,  Q  is  the  nominal  heat  flux  and  R  is  the  heat 
supply.  Lastly,  the  Clausius-Duhem  form  for  the  second  law  of  thermodynamics  is  expressed 
as 


P  :  F  +  Tt)  -  e  -  ^Grad  [T]-Q  .  (3.39) 

3.2  Linearized  Thermoelasticity 


We  consider  the  linearization  of  the  kinematic  quantities  from  the  previous  section;  for  the 
following  developments  see  Appendix  D  for  further  discussion.  Note  that  all  the  stress 
measures  are  equivalent  for  the  linearized  theory  and  we  will  therefore  use  a-  to  represent  the 
stress  tensor. 

Recall  the  local  form  of  the  balance  of  linear  momentum  may  be  expressed  as 


pi)  =  div  [crj  +  pb  in  Q  X  I 


(3.40) 


while  the  balance  of  angular  momentum  leads  to  symmetry  of  the  stress  tensor 


cr  =  in  X  /  (3.41) 

where  p  >  0  is  the  density,  v  is  the  acceleration  and  b  is  the  body  force.  The  local  form  of 
the  balance  of  energy  or  first  law  of  thermodynamics  may  be  expressed  as 


pe  =  —div  [g]  +  cr :  £  +  pr  in  Q  x  I  (3.42) 

where  e  is  the  internal  energy,  q  is  the  heat  flux,  e  is  the  infinitesimal  strain  tensor  and  r 
denotes  the  heat  supply.  The  local  form  of  the  second  law  of  thermodynamics  is 
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(3.43) 


pr]  +  div 


f 

-Pj;>0 


in  Q  X  I 


where  77  is  the  entropy  and  9  is  the  absolute  temperature.  The  Clausius-Duhem  inequality 
is  obtained  by  eliminating  the  heat  source  term,  r,  from  (3.42)  by  substituting  from  (3.43) 
to  obtain 


pTr]  —  pe  +  CT  :  e 

^  — "  V* 

'^loc 


-^grad  [T]  •  g  >  0 

' - V - ' 

T^con 


in  X  / 


(3.44) 


where  I?;oc  represents  the  local  or  internal  dissipation  and  the  dissipation  due  to  heat 
conduction  after  Truesdell  &  Noll  [1992].  Incorporating  a  generic  form  for  the  internal 
energy,  e  =  e(£,  77,  x)  we  obtain  the  time  derivative  of  the  internal  energy  as 


.  de  .  de 
e^  —  :e+—n+^ 


de 


dr]'  di 


(3.45) 


where  i  represents  a  general  set  of  internal  variables  representing  the  inelastic  response  of 
the  material.  Substituting  (3.45)  into  (3.44)  yields 


- 1)  ^ -  i)  -  -  i  • « -  I"!  ■  ’  ^ 

Assuming  Fourier’s  law  for  the  constitutive  equation  for  the  heat  flux,  q  =  —k  grad  [0],  the 
heat  conduction  dissipation  is  nonnegative,  Vcon  >  0.  Noting  that  (3.46)  is  linear  in  77,  e 
and  i,  independent  of  their  arguments  and  that  the  dissipation  inequality  must  hold  for  all 
admissible  processes  we  obtain 


de 

de 


de  . 

and  -  >  0. 

di 


(3.47) 
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We  introduce  an  alternative  form  for  the  balance  of  energy  for  subsequent  use  when  casting 
the  balance  equations  into  a  numerical  setting^.  Substituting  the  definition  for  Vioc  into 
(3.42)  results  in 


pTr]  = pr -\-Vioc  in  fix/.  (3.48) 

Introduce  the  Helmholtz  free  energy  density 


'*/’  =  '0(e^,  r,  a:)  =  e(e,  r),^]x)  -  Tr]  . 
Substituting  (3.49)  into  the  constitutive  equations  (3.47)  yields 


p7]=- 


dtp 


a  = 


dip 

de 


Taking  the  time  derivative  of  the  entropy  results  in 


(3.49) 


(3.50) 


f]  =  - 


dedT  ■  ^  dTdT  didT  '  ^ 


(3.51) 


Substituting  (3.51)  into  (3.48)  yields 


-PT =  -div  [q]  +  pr  +  Vioc  +  :  e  + 


dTdT~  ”■  ‘  de&T 

ct  =  -div  [q]  +  pr  +  Vtoc  -  H 


d$dT 


(3.52) 


^Typically,  algorithms  found  in  the  literature  for  coupled  thermomechanical  problems  consider  the  motion 
and  the  temperature  as  the  independent  variables. 
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where  c  =  >  0  is  referred  to  as  the  heat  capacity  and  H  =  T :  k  +  T  is  the 

structural  heating. 

Specializing  the  above  results  to  a  one-dimensional  bar,  neglecting  inertia  and  transient 
effects,  we  obtain 


o,x  +  pb  =  0  (3.53) 

-q,x  +  pr  +  T>ioc  -  /f  =  0  infix/ 

The  respective  boundary  conditions  for  (3.53)  are  given  as 

u  —  u  on  and  T  —  T  on  rr  (3.54) 

a  =  a  on  To-  and  q  —  q  on  F, 

where  we  have  decomposed  the  boundary  into  two  parts  for  each  problem 

r  =  r„ur,  =  Frur,  and  r„  n  =  Ft  n  r,  =  0  .  (3.55) 

When  supplemented  with  the  proper  boundary  conditions  (3.54),  the  balance  laws  (3.53) 
constitute  a  Boundary  Value  Problem. 


3.3  Constitutive  Equation 

To  complete  the  boundary  value  problem  for  the  linearized  theory  of  thermoelasticity  we 
shall  assume  that  the  one-dimensional  stress  response  function  has  the  form 


a  =  E[e-eL(e-  T)  -ot^T-  To)]  (3.56) 

where  a  is  the  stress,  E  is  the  elastic  modulus,  e  is  the  total  strain,  Sl  is  the  maximum 
residual  strain  obtained  by  detwinning  multiple  variant  martensite  (Bain  or  transformation 
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strain),  and  are  the  volume  fractions  of  the  positive  and  negative  variants  of  the 
martensite  twins  Achenbach  et.al.  [1986]  (see  Figure  3.1),  a.  is  the  coefiicient  of  thermal 
expansion  and  T  and  To  are  the  current  and  reference  temperature  fields  respectfully. 


Figure  3.1:  Variants  associated  with  various  lattice  configurations. 

Remark  3.1. 


1.  For  any  state,  the  simple  algebraic  relation  holds  +^a  =  Ij  where  represents 

the  volume  fraction  of  austenite  present.  Note  that  the  total  martensite  fraction  is  the 
sum  of  the  variants  ^  while  the  absence  of  both  and  indicate  the 

material  is  completely  austenite. 

2.  For  states  where  the  material  is  considered  to  be  in  a  self-accommodating  state 

in  which  that  material  is  100%  multiple  variant  martensite.  Although  its  configuration 
is  similar  to  100%  austenite  with  regard  to  the  overall  deformation  on  the  crystal  being 
zero,  the  crystal  structure  is  not  the  same. 

3.  For  a  one-dimensional  state  we  may  accurately  capture  the  behavior  of  the  phase 
transformation  using  two  internal  variables,  namely  and  For  higher  dimensions, 
the  underlining  physics  is  more  complicated  and  thus  additional  internal  variables  or 
variants  should  be  accounted  for;  see  Boyd  &  Lagoudas  [1996A,B]. 

4.  Different  authors  Tanaka  [1986],  Liang  k  Rogers.  [1990],  Brinson  [1993]  and  Brinson 
ET.  AL.  [1993]  have  utilized  linear  mixture  rules  for  the  elastic  material  modulus  E  = 
Ea  +  ^{Em  —  Ea)  and  coefficient  of  thermal  expansion  a  =  «„  +  ^{oim  —  cca),  where 
Ea,  Ejn,  cta  and  am  are  the  elastic  moduli  and  thermal  coefficient  of  expansion  for 
pure  austenite  and  martensite.  The  mixture  rules  enhance  the  model  by  accounting 
for  different  material  properties  at  each  pure  state.  The  introduction  of  these  mixture 
rules  may  or  may  not  be  justifiable  depending  upon  the  particular  alloy.  However,  the 
addition  of  these  mixtures  rules  presents  no  difficultly  in  the  present  formulation;  for 
simplicity  we  will  assume  the  moduli  to  be  constant. 
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3.4  Evolution  Equations 


In  the  following  subsections,  we  develop  explicit  expression  for  the  evolution  equations  for 
various  regions  within  the  phase  diagram. 

To  account  for  both  tensile  and  compressive  behavior,  we  consider  a  mapping  of  the  standard 
phase  diagram  in  the  tensile  stress  region  into  the  compressive  region.  For  the  present  work, 
the  transformation  parameters  are  denoted  with  either  a  or  a  “  representing  the  parameters 
in  the  tensile  and  compressive  regions  of  stress,  respectively.  The  ability  to  differentiate 
between  the  various  transformation  parameters  enables  the  model  to  qualitatively  capture 
the  behavior  observed  in  experimental  data. 


3.4.1  Production  of  Austenite 

Since  austenite  has  only  one  form  (variant),  it  is  sufficient  to  consider  the  evolution  of 
the  total  martensite  fraction.  The  positive  and  negative  variants  are  assumed  to  evolve 
proportional  to  the  total  martensite  fraction.  The  evolution  of  the  total  martensite  fraction 
may  be  expressed  in  an  integrated  form  as  a  linear  interpolation  between  the  start  and  finish 
transformation  lines  within  the  phase  transformation  region  shown  in  Figure  3.2 


Figure  3.2:  Admissible  transformation  regions  for  the  production  of  austenite.  The  points 
‘c’  and  ‘d’  denote  the  start  and  finish  temperatures  for  the  transformation,  while  Ks  and  Vaf 
are  the  start  and  finish  lines  for  the  transformation  and  the  directionality  of  the  evolution  is 
indicated  by  the  arrows. 

The  evolution  equation  is  expressed  as  follows: 
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The  parameter  is  introduced  to  account  for  cyclic  behavior  within  the  transformation  zone 
and  represents  the  maximum  value  of  ^  for  any  previous  loading  history.  The  parameters  Vas 
and  Vaf  denote  the  critical  values  of  stress  between  which  the  phase  transformation  occurs 
at  a  fixed  temperature,  where  the  subscripts  as  and  af  designate  stress  values  on  modified 
(by  the  initial  conditions)  austenite  start  and  finish  lines,  respectively.  The  starting  value  of 
the  phase  transformation  is  taken  to  be  a  function  of  the  initial  fractions  present,  while  the 
finish  value  is  taken  to  be  constant;  thus 

Vas  =  eras  +  Y - ^ 

t  -  qo 

Vaf  —  eTaf 

where  represents  the  initial  fraction  of  martensite  for  the  first  occurrence  of  the  trans¬ 
formation  and  aas  and  Ua/  are  defined  from  the  virgin  phase  transformation  lines  in  the 
stress-temperature  space  as 


aas  =  CaiT-Tas)  and  (Jaf  =  Ca{T-Taf)  (3.59) 

and  Ca  is  the  slope  of  the  transformation  lines,  assumed  fixed.  The  evolution  of  the  positive 
and  negative  variants  are  assumed  to  occur  in  proportion  to  their  existence  at  the  beginning 
of  the  phase  transformation 

e  =  and  r  =  ^e-  (3.60) 

sp  Cp 

Remark  3.2. 

•  To  accommodate  experimental  data  the  transformation  parameters  {Ca,  Ks,  Vaf)  may 
be  altered  to  reflect  differing  values  in  tension  and  compression. 
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State  and  Internal  Variables 


To  complete  the  model  for  the  production  of  austenite,  explicit  expressions  for  the  state 
and  internal  variables  are  derived  below.  In  addition,  material  tangents  will  also  be  derived 
in  order  to  facilitate  coding  the  model  into  a  numerical  setting,  such  as  the  finite  element 
method  where  a  tangent  matrix  is  required  for  a  solution  of  the  global  variables. 

If  the  phase  transformation  is  occurring,  the  inelastic  state  and  internal  variables  must  be 
determined.  Utilizing  the  constitutive  equation  (3.56)  and  the  evolution  equation  for  the 
production  of  austenite  (3.57),  we  may  solve  these  two  linear  equations  for  the  stress  and 
martensite  fractions  during  the  inelastic  process: 

a  =  |(£  -  £>■  -  -  {,-]  +  |)>„  -  1)^  -  c[T  -  T(,])  (3.61) 

where 

^  =  (3.62) 

whereas,  during  an  elastic  process  we  simply  have 

a  =  E[e-e^-  -  C]  -  o[T  -  To])  (3.63) 

To  obtain  the  mechanical  and  thermal  material  moduli  for  the  inelastic  case  we  take  the 
variation  of  the  stress  response  (3.61)i.  If  5  is  taken  as  the  variation  operator,  then: 

=  (i,„-l)^-Q  6T.  (3.64) 
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For  the  case  of  an  elastic  process  we  simply  obtain 


5a  =  —6s  +  —ST  =  E5e  -  aE5T  .  (3.65) 

3.4.2  Production  of  Positive  Single  Variant  Martensite 

For  the  production  of  positive  single  variant  martensite  it  is  sufficient  to  consider  the  evo¬ 
lution  of  one  of  the  variants,  namely  the  positive  variant.  The  remaining  negative  variant 
is  then  evolved  proportional  to  the  positive  variant  martensite.  The  evolution  of  the  sin¬ 
gle  positive  variant  martensite  fraction  may  be  expressed  in  an  integrated  form  as  a  linear 
interpolation  within  the  phase  transformation  shown  in  Figure  3.3: 


Figure  3.3:  Admissible  transformation  regions  for  the  production  of  positive  single  variant 
martensite.  The  point  ‘b’  denotes  the  temperature  below  which  austenite  is  not  stable. 
The  parameters  and  are  the  start  and  finish  lines  for  the  transformation  and  the 
directionality  of  the  transformation  is  indicated  by  the  arrow. 


{+  =  !  +  (!-  f+) 


-  KT/ 
KT/  -  k:. 


(3.66) 


The  parameter  is  introduced  to  account  for  cyclic  behavior  within  the  transformation  zone 
and  represents  the  maximum  value  of  for  the  previous  loading  history.  The  parameters 
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and  denote  the  critical  values  of  stress  between  which  the  phase  transformation 
occurs,  where  the  subscripts  ms  and  mf  designate  stress  values  on  a  modified  (by  the  initial 
conditions)  martensite  start  and  finish  lines,  respectively.  The  starting  value  of  the  phase 
transformation  is  taken  to  be  a  function  of  the  initial  fractions  present,  while  the  finish  value 
is  taken  to  be  constant;  thus 


v*.  =  (V*,  -  V*,)  (3.67) 

~  SO 

Vrnf  ~ 

yms  =  <^rns  +  +  [1  “  “  niin[^^,  ^o”])  (<^m/  “  Oms) 

where  represents  the  initial  fraction  of  martensite  for  the  first  occurrence  of  the  transfor¬ 
mation,  is  the  step  function  defined  as 


(7  <  0 
(7  >  0 


(3.68) 


and  Gms  and  Omf  are  defined  from  the  virgin  phase  transformation  lines  in  the  stress- 
temperature  space  as 


=  cFlr  +  C'm(T  -  Tms)  and  amf  =  olr  +  C'm(T  -  T^s)  (3.69) 


where  the  Macauley  bracket  (•)  is  defined  by 


a  <  0 
a  >  0 


(3.70) 


and  Cm  is  the  slope  of  the  transformation  lines,  assumed  fixed.  The  evolution  of  the  negative 
variant  is  assumed  to  occur  in  proportion  to  its  existence  at  the  beginning  of  the  phase 
transformation: 
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r  =  (3.71) 

Remark  3.3. 

•  In  equation  (3.67)3  the  second  term  on  the  right  hand  side  is  introduced  to  ensure  that 
only  the  single  variant  states  effect  the  initial  critical  stress. 

State  and  Internal  Variables 

The  model  is  completed  by  developing  expressions  for  the  state  and  internal  variables  and 
material  tangents  outlined  below. 

During  this  phase  transformation,  the  inelastic  state  and  internal  variables  must  be  deter¬ 
mined.  Utilizing  the  constitutive  equation  (3.56)  and  the  evolution  equation  for  the  produc¬ 
tion  of  positive  single  variant  martensite  (3.66),  we  may  explicitly  determine  the  stress  and 
martensite  fractions  during  the  inelastic  process  as 


6o  —  1  + 


V  KT/  -vz)  ■ 


whereas,  for  the  elastic  state  we  simply  have 


(3.73) 
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(3.74) 


<T  =  B  (e  -  e'  -  ei(5,+  -  f;|  -  a[T  -  T„]) 

5+  =  it  ■ 

To  obtain  the  mechanical  and  thermal  material  moduli  for  the  inelastic  case  we  take  the 
variation  of  the  stress  response  (3.72)i 


da  da  E  E 


ST 


(3.75) 


For  the  case  of  an  elastic  state  we  simply  obtain 


Sa  =  ^Se  +  ^5T  =  ESe  -  aEST  . 
de  dT 


(3.76) 


3.4.3  Production  of  Negative  Single  Variant  Martensite 

For  the  production  of  negative  single  variant  martensite  it  is  sufficient  to  consider  the  evo¬ 
lution  of  one  of  the  variants,  namely  the  negative  variant.  The  remaining  positive  variant 
is  then  evolved  proportional  to  the  single  negative  variant  martensite.  The  evolution  of  the 
single  negative  variant  martensite  fraction  may  be  expressed  in  an  integrated  form  as  a  linear 
interpolation  within  the  phase  transformation  shown  in  Figure  3.4. 

The  evolution  equation  is  expressed  as  follows 


IT  -  V, 


mf 


V~r  -  V- 

mf  ^  ms 


(3.77) 


The  parameter  is  introduced  to  account  for  cyclic  behavior  within  the  transformation  zone 
and  represents  the  maximum  value  of  for  the  previous  loading  history.  The  parameters 
and  V~j:  are  defined  as 
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Figure  3.4:  Admissible  transformation  regions  for  the  production  of  positive  single  variant 
martensite.  The  point  ‘b’  denotes  the  temperature  below  which  austenite  is  not  stable. 
The  parameters  and  are  the  start  and  finish  lines  for  the  transformation  and  the 
directionality  of  the  transformation  is  indicated  by  the  arrow. 


CT 


Figure  3.5:  Admissible  transformation  regions  for  the  production  of  negative  single  variant 
martensite.  The  parameters  and  V~j  indicate  the  start  and  finish  lines  for  the  transfor¬ 
mation  and  the  direction  for  the  production  of  single  variant  martensite  is  depicted  by  the 
arrow. 


V-,  =  (V-  _  V-,)  (3.78) 

SO 

Kn/  ~  ^rnf 

^ms  =  +  [1  -  n{a)](o  -  min[^^,^0-])  {amf  -  (Jms)  ■ 
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The  evolution  of  the  positive  variant  is  assumed  to  occur  in  proportion  to  its  existence  at 
the  beginning  of  the  phase  transformation 


i-r 


(3.79) 


State  and  Internal  Variables 

The  state  and  internal  variables  and  material  tangents  are  developed  below  to  complete  the 
model. 

During  this  phase  transformation,  the  inelastic  state  and  internal  variables  must  be  deter¬ 
mined.  Utilizing  the  constitutive  equation  (3.56)  and  the  evolution  equation  for  the  produc¬ 
tion  of  negative  single  variant  martensite  (3.77),  we  may  explicitly  determine  the  stress  and 
martensite  fractions  during  the  inelastic  process  as 


a  —  —  {  £  —  +  El  +  {bo  —  1) 


bo 


r  =  i  +  (i-e;) 


E 


a[T  -  To] 


mf 


-  y- 

^mf  ^  ms 


(3.80) 


where 


5o  =  1  ~ 


\v-f-v-, ) 


(3.81) 


whereas,  for  the  elastic  state  we  simply  have 


a  =  E{s-e^-  ex,[e;  -  ^^1  -  oc[T  -  To])  (3.82) 

C  =  Cp  ■ 
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To  obtain  the  mechanical  and  thermal  material  moduli  for  the  inelastic  case  we  take  the 
variation  of  the  stress  response  (3.80)i 


do  do  E  E 

5o  =  —6e  +  -^ST  =  —5e  +  — 
de  dT  bo  bo 


ibo-l)^{T-Tms)-a 


5T 


(3.83) 


3.4.4  Production  of  Multiple  Variant  Martensite 


Since  multiple  variant  martensite  has  equally  distributed  proportions  of  positive  and  negative 
variants  of  martensite,  from  a  zero  initial  state"^,  it  is  sufficient  to  consider  the  evolution  of 
the  total  martensite  fraction.  The  total  positive  and  negative  variants  are  then  evolved 
proportional  to  the  total  martensite  fraction  present.  The  evolution  of  the  total  martensite 
fraction  may  be  expressed  in  an  integrated  form  as  a  linear  interpolation  within  the  phase 
transformation  zone  shown  in  Figure  3.6. 


CT 


Figure  3.6:  Admissible  transformation  regions  for  the  production  of  multiple  variant  marten¬ 
site.  The  points  ‘a’  and  ‘b’  denote  the  start  and  finish  temperatures  for  the  phase  transfor¬ 
mation  and  the  direction  of  the  evolution  is  depicted  by  the  arrow. 

The  evolution  equation  for  the  production  of  multiple  variant  martensite  is  expressed  as 

«  =  l  +  ■  (3.84) 

\”m/  ”ms  / 

'^Producing  multiple  variant  martensite  from  100  %  austenite. 
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The  parameter  is  introduced  to  account  for  cyclic  behavior  within  the  transformation  zone 
and  represents  the  maximum  value  of  (  for  the  previous  loading  history.  The  parameters  0ms 
and  Omf  denote  the  critical  values  of  temperature  between  which  the  phase  transformation 
occurs,  where  the  subscripts  ms  and  mf  designate  temperature  values  on  modified  (by  the 
initial  conditions)  martensite  start  and  finish  lines,  respectively.  The  starting  value  of  the 
phase  transformation  is  taken  to  be  a  function  of  the  initial  fractions  present,  while  the  finish 
value  is  taken  to  be  constant;  thus 


0ms  =  Tms  +  (Tmf  "  Tms)  (3.85) 

Omf  ~  Tmf 

where  represents  the  initial  fraction  of  martensite  for  the  first  occurrence  of  the  transfor¬ 
mation  and  Tms  ^md  Tmf  are  material  parameters  taken  from  the  virgin  phase  transformation 
lines  in  the  stress-temperature  space. 

The  evolution  of  the  positive  and  negative  variants  are  assumed  to  occur  in  proportion  to 
their  existence  at  the  beginning  of  the  phase  transformation 


=  +  and  r  =  f;  +  5(«-fp)- 


(3.86) 


State  and  Internal  Variables 

The  model  is  completed  by  deriving  expressions  for  the  state  and  internal  variables  and 
material  tangents  outlined  below. 

During  this  phase  transformation,  the  inelastic  state  and  internal  variables  must  be  deter¬ 
mined.  Note  that  the  evolution  equation  for  the  production  of  multiple  variant  martensite 
(3.84)  is  independent  of  the  stress  field,  hence  we  may  simply  calculate  the  martensite  fraction 
(assuming  a  known  temperature)  and  substitute  those  results  into  the  constitutive  equation 
(3.56)  as 
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(3.87) 


^  =  1  +  (1  - 


T  —  9m  f  \ 

^mf  ^ms  / 


a  =  E  {e  -  -  SLlt  -  n  -  €^[T  -  To])  . 


To  obtain  the  mechanical  and  thermal  material  moduli  for  the  inelastic  case  we  take  the 
variation  of  the  stress  response  (3.87)2 


5a  =  ^5e  +  ^5T  =  E5e  -  aE5T  . 
oe  oT 


(3.88) 
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Chapter  4 

Finite  Element  Developments 


This  section  covers  the  algorithmic  approximation  of  the  constitutive  theory  in  the  setting 
of  the  finite  element  method.  Due  to  the  nature  of  the  nested  elastic  and  inelastic  regions, 
the  determination  of  when  a  transformation  is  active  becomes  paramount  to  the  model’s 
implementation.  This  detection  procedure  is  known  as  state  determination.  An  outline  of 
the  algorithmic  procedure  is  described  and  discussed  below. 

Lastly,  to  examine  the  properties  of  the  resulting  equations  the  model  is  implemented  within 
the  context  of  line  elements  in  the  finite  element  method.  Specifically,  we  numerically  assess 
the  behavior  of  a  class  of  shape  memory  alloys  with  linear  and  non-linear  kinematics  for 
multi-dimensional  truss  bars  and  a  two-dimensional  beam. 


4.1  State  Determination 


In  this  section  we  outline  the  procedure  used  for  state  determination  at  time  tn+i  using 
information  from  the  current  global  iterate  and  the  previous  time  step  (•)„  for  various 

phase  transformations.  Also  note  we  consider  that  the  temperature  is  given  at  time  by 
either  a  prescribed  load  history  or  solving  the  heat  conduction  problem. 

The  phase  space  diagram  in  Figure  2.4  has  lead  many  researcher  to  use  existing  schemes 
for  state  determination,  such  as  trial  state  methods  in  classical  plasticity.  Unfortunately, 
classical  algorithms  fail  to  determine  the  correct  state  even  for  simple  simulations  due  to 
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the  wide  range  of  trial  states.  As  an  example,  during  the  loading  portion  of  a  pseudoelastic 
process,  both  elastic  and  inelastic  states  are  admissible,  since  the  inelastic  region  is  bounded 
by  two  elastic  regions.  Simply  checking  the  trial  value  of  stress  with  the  “yield  function” 
and  its  directionality  will  not  guarantee  the  correct  state. 

We  introduce  a  new  algorithm  which  accounts  for  nested  elastic  and  inelastic  regimes.  The 
underlining  concept  relies  on  the  use  of  modified  trial  state  variables  and  the  determination 
of  trial  states  within  an  admissible  transformation  region.  The  resulting  trial  states  are 
checked  for  consistency  with  respect  to  the  phase  diagram,  i.e.  directionality  and  magnitude 
of  the  stress,  temperature  and  internal  variables  (martensite  fractions). 


4.1.1  Production  of  Austenite 

Consider  the  production  of  austenite  with  initial  conditions  at  time  shown  in  Figure  4.1 
Given  a  temperature  increase  to  r„+i  we  wish  to  determine  the  state  for  time 


Figure  4.1:  Bounds  for  martensite  fraction  in  phase  space.  From  we  may  develop  an 
isofraction  line  on  the  phase  diagram,  below  which  transformations  are  elastic. 

For  the  production  of  austenite,  we  consider  the  state  (T„,a„)  on  the  phase  diagram  in 
Figure  4.2  and  two  critical  states  (Tn+i,(7i)  and  {Tn+i,a2)  at  time  tn+i-  The  pair  (r„+i,cri) 
corresponds  to  a  point  on  the  iso-fraction  line,  while  the  pair  (T„+i,cr2)  corresponds  to  the 
point  at  which  the  martensite  transformation  is  completed.  From  these  two  points  we  see  at 
time  tn+i  there  exist  three  possible  regions  for  the  solution,  see  Figure  4.2. 

To  determine  the  correct  trajectory  for  the  state  at  we  evaluate  the  constitution  for  all 
three  regions.  The  three  possible  choices  are: 


2-45 


Figure  4.2:  Possible  trajectories  in  stress-fraction  space.  We  see  that  values  of  stress  greater 
than  (72  and  less  than  cti  result  in  elastic  states. 


1.  (cr  <  (7i  -  Elastic  State)  Evaluate  the  constitution  using  the  previous  converged  value 
for  the  internal  variables: 

=  £  (4+V’  -  <  -  IC  -  C]  -  [r»+i  -  Hi)  (4.1) 

2.  ((7i  <  (7  <  (72  -  Production  of  Austenite)  Simultaneously  solve  for  the  constitution  and 
the  evolution: 

<7=1  (4+V’  -  <  -  -  (;j  +  [6.  -  IJ^  -  alH«  -  T„l)  (4.2) 

f  =  £  -f  ( 


3.  (a  >  (72  -  Elastic  State)  Evaluate  the  constitution  assuming  that  the  martensite  trans¬ 
formation  is  completed: 


where 


■^  =  £  [C  -  Cp.]  -  a  IH+i 


(4.3) 


{&- 


^apx  —  max 


^apx  ^p  /  ^p^dpx 
^apx  /  ^p^apx 


^0 


{T-Tas) 
Taf  -  Tas 


(4.4) 


Remark  4.1. 


•  The  fraction  ^apx  represents  the  maximum  amount  of  martensite  which  may  be  trans¬ 
formed  into  austenite  for  a  particular  temperature  T  and  initial  conditions 

For  the  production  of  austenite  to  occur  T  >Tas  +  {l  —  ^o){Taf  —  Tas)  and  >  ^apx-  If  this 
criteria  is  met  then  paths  1,  2  and  3  are  evaluated.  The  additional  consistency  checks  on 
the  directionality  and  magnitudes  of  the  stress  and  the  internal  variables  for  each  path  are 
outlined  below. 


Path  1 

•  If  <7  >  0  and  <7  <  cTas  +  “  <7as)  then  admissible. 

•  If  (T  <  0  and  <7  <  —  (7a5  —  (1  —  ^n)((7af  —  <7 as)  then  admissible. 

Path  2 

•  If  (7  >  0,  cr  <  I/+,  —da+C^ dT  >  0  and  ^  then  admissible,  where  da  =  a— an 
and  dT  =  Tn-\-i  —  T„. 

•  If  a  <  0,  V~g  <  a,  da  —  C~dT  >  0  and  ^  then  admissible. 

Path  3 

•  If  (7„  >  0  and  a  <  max[0,  then  admissible. 

•  If  (Tn  <  0  and  a  >  min[0,  V~j]  then  admissible. 


If  multiple  states  are  admissible  we  select  the  state  which  is  closest  to  the  previous  converged 

state  via  distance  metric  in  stress  space.  By  choosing  the  closest  state  to  the  previous  state 

we  restrict  the  local  behavior  of  the  local  constitution  from  returning  spurious  states.  The 

metric  is  minimized  using  the  admissible  current  states  and  the  previous  converged  state,  i.e. 

cr  =  min[d((7„,  cTj)]  where  d{x,y)  —\x  —  y\.  The  resulting  state  is  returned  to  the  element 
i 

and  the  residual  and  tangent  arrays  are  built. 
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Figure  4.3:  Bounds  for  martensite  fraction  in  phase  space.  The  parameters  corresponds 
to  the  minimum  value  of  stress  required  to  initiate  a  transformation  and  the  arrow  indicates 
the  direction  of  the  transformation. 

4.1.2  Production  of  Positive  Single  Variant  Martensite 

Consider  now  the  production  of  positive  single  variant  martensite  with  initial  conditions  at 
time  tn  shown  in  Figure  4.3.  Given  a  temperature  increase  to  Tn+i  we  wish  to  determine  the 
state  for  time 

For  the  production  of  single  variant  positive  martensite  we  consider  the  state  (T„,  cr„)  on  the 
phase  diagram  in  Figure  4.3  and  two  critical  states  {Tn+\,ai)  and  (Tn+i,a2)  at  time  tn+i- 
The  pair  corresponds  to  a  point  on  the  iso-fraction  line,  while  the  pair  (T„+i,cr2) 

corresponds  to  the  point  at  which  the  martensite  transformation  is  completed.  From  these 
two  points  we  see  at  time  tn+i  there  exists  three  possible  regions  for  the  solution,  see  Figure 
4.4. 

To  determine  the  correct  trajectory  for  the  state  at  tn+i  we  evaluate  the  constitution  for  all 
three  regions.  The  three  possible  choices  are: 

1.  {a  <  ai  -  Elastic  State)  Evaluate  the  constitution  using  the  previous  converged  value 
for  the  internal  variables: 

<^  =  E  (4+Y’  K  -  c]  -  a  [T„+1  -  r„l)  (4.5) 
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Figure  4.4:  Possible  trajectories  in  stress-fraction  space.  We  observe  that  values  of  stress 
greater  than  <72  and  less  than  (j„  result  in  elastic  states. 


2.  <  o  <  02  -  Production  of  Single  Variant  Martensite)  Simultaneously  solve  for  the 

constitution  and  the  evolution: 


E  I  _(fc+i) 


a  =  -  1 4Tr  -  -eL  +  {bo  -  1)^  -  -  To] 


(4.6) 


r  =  i+(i-^;j  ^ 


^  -  v^f 


y+.  _  7+ 

^  mf  ^  ms  ^ 

3.  3B  (a  >  (72  -  Elastic  State)  Evaluate  the  constitution  assuming  that  the  martensite 
transformation  is  completed: 


<’  =  E  (4+V’  -si- CL -a  17;+,  -  Til)  (4.7) 

For  the  production  of  positive  single  variant  martensite  to  occur  <  1  and  (7„  >  0.  If  this 
criteria  is  met  then  paths  1,  2  and  3  are  evaluated.  The  additional  consistency  checks  on 
the  directionality  and  magnitudes  of  the  stress  and  the  internal  variables  for  each  path  are 
outlined  below. 


Path  1 


•  If  O’  <  ams  +  ^i{omf  —  Oms)  then  admissible. 


Path  2 


•  If  (7  >  Oms  +  (^m/  ~  Oms),  do  —  C^dT  >  0  and  then  admissible. 
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Path  3 


•  If  cr  >  then  admissible. 


If  multiple  states  are  admissible  we  select  the  state  which  is  closest  to  the  previous  converged 
state  via  the  metric  a  =  min[d((Tn,  aj)].  The  resulting  state  is  returned  to  the  element  and 

i 

the  residual  and  tangent  arrays  are  built. 


4.1.3  Production  of  Negative  Single  Variant  Martensite 

Consider  now  the  production  of  negative  single  variant  martensite  with  initial  conditions  at 
time  tn  shown  in  Figure  4.5.  Given  a  temperature  increase  to  we  wish  to  determine  the 
state  for  time  . 


Figure  4.5:  Bounds  for  martensite  fraction  in  phase  space.  Note  for  the  production  of  single 
variant  martensite  to  occur  the  stress  must  be  greater  than  cr„  and  less  than  02  and  the 
resulting  state  must  be  evolving  in  the  direction  of  the  arrow. 

For  the  production  of  single  variant  negative  martensite  we  consider  the  state  (T„,  an)  on  the 
phase  diagram  in  Figure  4.5  and  two  critical  states  (T„+i,cri)  and  (T„+i,a2)  at  time  tn+i- 
The  pair  (Tn+i,cri)  corresponds  to  a  point  on  the  iso-fraction  line,  while  the  pair  {Tn+i,a2) 
corresponds  to  the  point  at  which  the  martensite  transformation  is  completed.  From  these 
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Figure  4.6:  Possible  trajectories  in  stress-fraction  space.  We  observe  from  the  diagram  that 
elastic  states  exist  for  a  >  an  or  a  <  a2.  For  values  of  stress  within  these  elastic  state  exist 
the  possibility  of  a  phase  transformation. 

two  points  we  see  at  time  tn+i  there  exists  three  possible  regions  for  the  solution,  see  Figure 
4.6. 

To  determine  the  correct  trajectory  for  the  state  at  tn+i  we  evaluate  the  constitution  for  all 
three  regions.  The  three  possible  choices  are; 


1.  (a  >  ai  -  Elastic  State)  Evaluate  the  constitution  using  the  previous  converged  value 
for  the  internal  variables: 

=  -B  (4+V’  -  eS  -  £1  [fi'  -  C]  -  o  lr„+i  -  Tol)  (4.8) 


2.  (ai  >  a  >  a2  -  Production  of  Single  Variant  Martensite)  Simultaneously  solve  for  the 
constitution  and  the  evolution: 

"  =  I  -e^+si.  +  (bo-l)^- oir„+,  -  To])  (4.9) 


\  ^  ms 


3.  ((7  <  (72  -  Elastic  State)  Evaluate  the  constitution  assuming  that  the  martensite  trans¬ 
formation  is  completed: 


<^  =  B  (4+?*  -el-n-a  [T„+,  -  T„]) 


(4.10) 
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For  the  production  of  negative  single  variant  martensite  to  occur  <  1  and  a„  <  0.  If  this 
criteria  is  met  then  paths  1,  2  and  3  are  evaluated.  The  additional  consistency  checks  on 
the  directionality  and  magnitudes  of  the  stress  and  the  internal  variables  for  each  path  are 
outlined  below. 

Path  1 


•  If  O’  >  -ams  -  “  <^ms)  then  admissible. 

Path  2 

•  If  cr  <  -ams  -  -  (^ms),  do  -  C~dT  <0  and  then  admissible. 


Path  3 


•  If  cr  <  then  admissible. 


If  multiple  states  are  admissible  we  select  the  state  which  is  closest  to  the  previous  converged 
state  via  the  metric  a  =  min[(f(an,  cTj)].  The  resulting  state  is  returned  to  the  element  and 

i 

the  residual  and  tangent  arrays  are  built. 


4.1.4  Production  of  Multiple  Variant  Martensite 

Consider  now  the  production  of  multiple  variant  martensite  with  initial  conditions  at  time 
tn  shown  in  Figure  4.7.  Given  a  temperature  increase  to  T„+i  we  wish  to  determine  the  state 
for  time  tn+i- 

For  the  production  of  multiple  variant  martensite  we  consider  the  state  on  the 

phase  diagram  in  Figure  4.7  and  two  critical  states  (T„+i,cri)  and  (T„+i,(T2)  at  time  tn+i- 
The  pairs  (r„+i,cri)  and  {Tn+i,a2)  corresponds  to  limit  points  between  which  a  multiple 
variant  martensite  transformation  occurs.  Prom  these  two  points  we  see  at  time  there 
exists  one  possible  region  for  the  solution,  see  Figure  4.8,  provided  dT  <  0. 

To  determine  the  correct  trajectory  for  the  state  at  tn+i  we  evaluate  the  constitution  and 
evolution  for  this  region: 
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Figure  4.7:  Bounds  for  martensite  fraction  in  phase  space.  For  the  production  of  multiple 
variant  martensite  to  occur  >  o  >  ,  and  T  <  0. 
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Figure  4.8:  Possible  trajectories  in  stress-fraction  space.  Note  that  during  this  transforma¬ 
tion  the  evolution  is  decoupled  from  the  constitution  and  hence  only  an  evaluation  is  required 
to  determine  the  state  and  internal  variables. 


•  (<^2  >  c’’  >  cn  -  Production  of  Multiple  Variant  Martensite)  Simultaneously  solve  for 
the  constitution  and  the  evolution: 


«=i+(i-4) 

‘^  =  E  (e'yy  -  e'  -  ulC  -  ?„1  -  alT„+i  -  To]) 


(4.11) 


The  resulting  state  is  returned  to  the  element  and  the  residual  and 


tangent  arrays  are  biiilt. 
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4.1.5  Outline  of  Algorithm 


The  previous  sections  outline  the  various  constitutive  and  evolution  equations  which  need 
to  be  evaluated.  To  minimize  the  computational  effort,  we  may  exclude  states  by  simple 
conditional  checks  with  the  previous  state.  Once  a  reduced  set  of  states  is  determined, 
additional  checks  are  performed  to  ensure  that  the  directionality  and  magnitude  of  the 
resulting  state  is  consistent  with  its  location  on  the  phase  diagram. 

We  now  consider  the  various  transformations  which  can  occur  and  their  associated  consis¬ 
tency  check. 


Production  of  Austenite 

For  the  production  of  austenite  to  occur  T  >  Ta^  -f  (1  —  ^o)(T’a/  —  and  >  ^apx  must  be 
met.  If  this  criteria  is  met  then  paths  1,  2  and  3  are  evaluated.  The  additional  consistency 
checks  on  the  directionality  and  magnitudes  of  the  stress  and  the  internal  variables  for  each 
path  are  outlined  below. 

Path  1:  Elastic 

•  If  cr  >  0  and  a  >  then  admissible. 

•  If  (7  <  0  and  a  <  V~  then  admissible. 

Path  2:  Transformation 

•  If  <7  >  0,  a  <  -da  +  C^dT  >  0  and  ^apx  <  ^  <  then  admissible. 

•  If  <j  <  0,  V~  <a,  da  -  C~dT  >  0  and  ^apx  <^<in  then  admissible. 

Path  3:  Transformation  with  overshoot 

•  If  cr„  >  0  and  a  <  max[0,  then  admissible. 

•  If  (T„  <  0  and  a  >  min[0,  V~f]  then  admissible. 


2-54 


Production  of  Positive  Single  Variant  Martensite 


For  the  production  of  positive  single  variant  martensite  to  occur  <  1  and  cr„  >  0  must  be 
met.  If  this  criteria  is  met  then  paths  1,  2  and  3  are  evaluated.  The  additional  consistency 
checks  on  the  directionality  and  magnitudes  of  the  stress  and  the  internal  variables  for  each 
path  are  outlined  below. 


Path  1:  Elastic 

•  If  cr  <  then  admissible. 

Path  2:  Transformation 

•  If  a  >  da  —  dT  >  0  and  1  >  ^'''  >  ^^  then  admissible. 

Path  3:  Transformation  with  overshoot 

•  If  <7  >  then  admissible. 


Production  of  Negative  Single  Variant  Martensite 

For  the  production  of  negative  single  variant  martensite  to  occur  <  1  and  <  0  must  be 
met.  If  this  criteria  is  met  then  paths  1,  2  and  3  are  evaluated.  The  additional  consistency 
checks  on  the  directionality  and  magnitudes  of  the  stress  and  the  internal  variables  for  each 
path  are  outlined  below. 

Path  1:  Elastic 

•  If  cr  >  then  admissible. 

Path  2:  Transformation 

•  If  cr  <  da  —  C~dT  <  0  and  1  >  then  admissible. 

Path  3:  Transformation  with  overshoot 

•  If  cr  <  V~^  then  admissible. 
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Production  of  Multiple  Variant  Martensite 


For  the  production  of  multiple  variant  martensite  to  occur  <  1  must  be  met.  If  this 
criteria  is  met  then  the  path  ID  is  evaluated.  The  additional  consistency  checks  on  the 
directionality  and  magnitude  of  the  stress  temperature  and  the  internal  variables  for  the 
path  is  outlined  below. 


Path  2:  Transformation 

•  If  V~g  <  a  <  dT  <  0,  and  1  >  ^  >  then  admissible. 


Multiple  States 

Once  the  respective  constitution  and  evolution  equations  have  been  evaluated  and  consis¬ 
tency  checks  performed  the  state  is  returned  to  the  main  driver  element  for  the  construction 
of  the  residual  and  tangent  arrays  described  in  the  following  sections.  If  multiple  states 
are  admissible,  we  select  the  state  which  is  closest  to  the  previous  converged  state^  and  is 
implemented  with  a  metric.  The  metric  is  minimized  using  the  admissible  current  states  and 
the  previous  converged  state,  i.e.  a  =  min[d(cr„, (7^)]  where  d{x,y)  =\  x  —  y  \.  The  resulting 

i 

state  is  returned  to  the  element  and  arrays  built. 


4.2  Multi- Dimensional  Truss-Bar  Element 


In  this  section  we  consider  the  formulation  of  a  multi-dimensional  bar  element  for  both 
linear  and  nonlinear  kinematics,  see  Figure  4.9  for  a  representation  of  the  geometry.  In 
addition  to  the  physically  nonlinear  mechanical  response  of  the  material,  we  also  consider 
the  thermomechanical  effects  for  the  material.  We  begin  with  a  brief  outline  of  the  basic 
notation  used  throughout  the  developments  and  then  proceed  with  a  description  of  the  strain 
measure  used.  A  variational  equation  for  elastostatics  is  reviewed  and  the  finite  element 
interpolations  for  a  line  element  are  presented.  Lastly,  the  finite  element  arrays  necessary 
for  implementation  are  developed. 

^By  choosing  the  closest  point,  we  restrict  the  local  behavior  of  the  local  constitution  from  returning 
spurious  states. 
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Figure  4.9:  Reference  configuration  for  the  bar  element. 

4.2.1  Notation 

We  consider  an  initially  straight  bar  of  length  L  and  cross  section  ^  with  smooth 
boundary  9r2o  to  represent  a  bounded  reference  configuration  B  for  the  continuum  body.  We 
admit  the  decomposition  of  the  boundary  into  two  parts:  r„  C  dB  where  the  displacement 
is  prescribed  as  u  =  u  and  Fj  C  dB  where  the  traction  is  prescribed  as  =  d  subject  to 


=  and  r„nrt  =  0  (4.12) 

where  cr^  is  the  stress  directed  along  the  axis  or  principal  direction  of  the  bar.  For  subsequent 
treatment  of  the  variational  formulation  we  distinguish  two  classes  of  functions,  namely,  the 
space  admissible  solutions  and  the  space  of  admissible  variations. 

Let  U  be  the  space  of  admissible  displacements  written  as 


U  —  {u  \  u  E  H^{L)  and  u  =  u  on  Fu}  (4-13) 

and  V  be  the  space  of  admissible  displacement  variations  written  as 


V={5u\SuE  H^{L)  and  (5u  =  0  on  r„}  (4.14) 

where  is  the  Sobolev  space  of  degree  1,  consisting  of  functions  which  posses  square- 
integrable  first  order  derivative  and  are  themselves  square-integrable. 
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4.2.2  Strain  Measure 


The  infinitesimal  increment  of  the  strain  due  to  a  change  in  the  length  is  expressed  as 


ds,  =  j  (4.15) 

where  deg  is  the  increment  in  the  strain  directed  along  the  principal  axis  of  the  bar,  dl  is 
the  incremental  change  in  length  of  the  bar  and  I  is  the  deformed  length  of  the  bar.  By 
integrating  the  incremental  strain  from  the  initial  length  L  to  the  deformed  length  I  we 
obtain  the  total  strain 


e,  =  de,  =  In  =  In  =  ln(l  +  e)  (4.16) 

where  AL  represents  the  change  in  length  of  the  bar.  The  linear  strain  measure  may  be 
obtained  by  a  series  expansion  of  (4.16),  retaining  only  the  linear  terms  of  the  expansion. 
Expanding  the  last  expression  in  (4.16)  results  in 


g2  g3  g4 

Eg  =  ln(l  —  - =  e  +  o(e^) 

subject  to  — 1  <  e  <  1.  From  (4.17)  the  linearized  strain  measure  becomes 


(4.17) 


AL  dUg 


(4.18) 


where  Ug  is  the  deformation  directed  along  the  principal  axis  of  the  member  and  s  is  a 
distance  measure  directed  along  the  principal  axis  of  the  member. 
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4.2.3  Variational  Formulation 


An  approximate  solution  of  the  boundary  value  problem  is  constructed  from  a  variational 
statement  of  the  problem.  The  basic  field  equations  may  be  included  in  a  variational  state¬ 
ment  for  the  elasticity  problem  using  a  potential  functional.  Accordingly,  we  have  a  func¬ 
tional  in  which  the  displacement  field,  Uj  G  W  is  regarded  as  the  independent  variable.  The 
proposed  functional  IT  :  Z//  — )•  R  may  be  expressed  as 

n(u,)  =  [  W{es)  dV  +  Ue,t{us)  (4.19) 

Jb 

where  ^(£5)  is  a  stored  energy  function  and  for  conservative  external  loading 

UeM  =  -  fb.udV-  [  au  dA  (4.20) 

Jb  Jvt 

bs  being  the  body  force  per  unit  volume. 

We  may  state  the  problem  as:  Find  Ug  eU  which  makes  the  functional  11  (uj  stationary  for 
all  admissible  variations  Sug  €  V. 

The  stationary  point  of  fl  is  obtained  by  setting  to  zero  the  first  variation  of  (4.19)  with 
respect  to  the  independent  field.  Accordingly, 

SU=  [  ^dsg  dV-  [  Subg  dV-  [  Sua  dA  =  0  in  B  (4.21) 

Jb  Jb  Jvt 

for  all  admissible  variations  5ug  G  V. 

We  may  recast  the  first  term  of  (4.21)  as 

=  J  ^^0  ds  J  degUg  CIq  ds  =  SegOg  Vt^L  (4.22) 

which  represents  the  internal  virtual  work. 


2-59 


4.2.4  Finite  Element  Interpolations 


To  enable  numerical  implementation  of  the  variational  formulation  a  discretization  with 
respect  to  the  global  variables  is  performed.  We  begin  by  discretizing  the  bar  into  a  finite 
number  of  points  called  nodes.  Connecting  two  nodes  constitutes  an  element  on  which  the 
formulation  is  developed.  In  the  following  the  interpolations  for  the  global  variables  are 
introduced. 

The  selection  of  the  displacement  interpolation  functions  for  the  element  must  satisfy  certain 
basic  requirements  such  that  the  solution  converges  to  the  exact  solution  under  certain  special 
cases. 

The  displacement  interpolation  functions  must 


1.  be  continuous  within  the  element, 

2.  be  able  to  capture  rigid  body  modes, 

3.  be  able  to  capture  constant  strain  states, 

4.  ensure  inter-element  compatibility. 


The  bar  element  used  for  the  subsequent  developments  consists  of  an  initially  straight  bar 
with  two  nodal  points  located  on  the  boundary  of  the  line  segment.  The  obvious  selection 
for  the  interpolation  functions  which  satisfy  the  above  criteria  are 


iVi  =  1  -  y  and  N2  =  y  .  (4.23) 

Ju  Lj 

Utilizing  the  interpolation  functions  in  (4.23)  we  may  express  the  displacement  vector  u,  as 
well  as  the  position  vector  X  for  the  bar  as 


2 

Uf(5)  = 

1=1 


2 

and  -^1(5)  =  Y^NjXn 
1=1 


(4.24) 
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where  z  =  1,  •  •  •  ,  ndm  is  the  dimension  of  the  space  and  u  and  X  are  the  nodal  quantities  of 
the  displacement  and  coordinates,  respectively.  For  further  elaboration  on  the  development 
of  finite  element  interpolation  see  Hughes  [1987]  or  ZIENKIEWICZ  &  Taylor  [1989]. 


4.2.5  Linear  Kinematic  Element 

We  proceed  by  defining  the  strain  measure  used  for  the  formulation,  develop  the  expression 
for  the  internal  virtual  work  from  which  the  residual  equations  may  be  expressed.  Lastly, 
the  resulting  nonlinear  residual  equation  is  linearized  to  aff’ord  a  solution  of  the  nonlinear 
set  of  equations  via  Newton’s  method. 


Strain  Measure 

Recall  from  the  previous  section  the  expression  for  the  linearized  strain  measure  is 


^5 


TldTTL 

Ullg  X  ^  j  OUi 
1=1 


(4.25) 


where  Ui  is  the  displacement  in  the  coordinate  direction  i  and  li  is  the  direction  cosine  in  the 
coordinate  direction  defined  in  terms  of  the  nodal  locations  as 


V  V  A  V  nam 

=  -^  and  ^  (X2i  -  Xuf  •  (4.26) 

1—1 

Utilizing  (4.26)  we  may  express  the  strain  in  the  principal  direction  as 


AXi  Xui 
L  L 


(4.27) 


where  summation  over  i  is  implied  and  Auj  =  U2i  —  uu  represents  the  relative  displacement 
in  the  coordinate  direction. 
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Residual  and  Tangent  Arrays 


Substitution  of  (4.27)  into  the  expression  for  the  internal  virtual  work  yields 


^^int  — 


AXi  5Aui 

~T  IT 


—  ^AVtiCXgli^Q 


(4.28) 


or  in  matrix  notation  we  have 


^^int  —  ■{ 


(4.29) 


where  Fi  =  agh^o  is  the  internal  force. 

To  introduce  one-way  thermal  effects^  into  the  formulation  we  introduce  an  additional  nodal 
quantity  or  degree  of  freedom,  namely  temperature  into  (4.29) 


Su2i 


f 

ST2i  }  < 


(4.30) 


The  vector  post-multiplying  the  virtual  displacements  and  temperatures  in  (4.30)  is  com¬ 
monly  referred  to  as  the  internal  load  vector  and  is  denoted  as  Fint  for  future  developments. 

Remark  4.2. 


Recall  from  Chapter  3  that  the  stress  ag  is  in  general  a  function  of  the  mechanical  and 
thermal  state  and  internal  variables,  i.e.  ag  =  a(e,  T,  ^).  Hence,  its  variation  with  respect  to 
the  added  degree  of  freedom  is  nonzero  in  general  and  must  be  considered  when  constructing 
the  solution. 

^The  one-way  coupling  used  throughout  refers  to  the  ability  of  a  thermal  gradient  to  produce  mechanical 
strains,  but  the  reverse  process  is  uncoupled  and  hence  not  considered. 


2-62 


To  facilitate  a  solution  of  the  mixed  boundary  value  problem  the  nonlinear  equation  (4.21) 
is  linearized  and  solved  by  a  Newton’s  method  as  a  sequence  of  linearized  problems.  Hence, 
linearizing  (4.21),  neglecting  external  work  terms  we  obtain 


d((5n)  =  6es  dos  fio  L  (4.31) 

where 


d(7g  —  Dj^dSg  +  D'ldT  (4.32) 

and  Dm  and  Dt  are  the  mechanical  and  thermal  material  moduli.  Substituting  (4.32)  into 
(4.31)  with  the  aid  of  (4.27)  yields 


ndm  ndm 


d((5n)  =  EE 


SAudi 


i=l  j=l 

and  expressed  in  matrix  notation  we  have 


^  dAujlj  ^  dTj 


L 


(4.33) 


<i(OT)  = 


where 


(  duu 

T 

kij 

hi 

—  k 

rvij 

hi 

STu 

0 

0 

0 

0 

f 

l5«2i 

> 

—kij 

hi 

k  •  ■ 

hi 

mi  J 

0 

0 

0 

0 

kij  - 

=  h 

T. 

and 

hi  — 

(  duij  'I 
dTij 
dU2j 
dT2j 


>  =  St]^  k  dr} 


(4.34) 


f2o 

‘T 


(4.35) 


The  matrix  post-multiplying  the  virtual  displacements  and  temperatures  in  (4.34)  is  com¬ 
monly  referred  to  as  the  tangent  matrix  and  is  denoted  as  k.  The  terms  Dm  and  Dt  are 
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the  tangent  moduli  from  the  previous  sections  depending  on  the  state  of  the  material.  Note 
that  the  continuum  and  algorithmic  tangents  are  the  same  since  the  evolution  of  the  internal 
variables  is  treated  exactly. 


4.2.6  Non-Linear  Kinematic  Element 

We  proceed  by  defining  the  strain  measure  used  for  the  formulation  and  then  develop  the 
expression  for  the  internal  virtual  work  from  which  the  residual  equations  may  be  expressed. 
Lastly,  the  resulting  nonlinear  residual  equation  is  linearized  to  allow  solution  of  the  nonlinear 
set  of  equations  via  Newton’s  method. 


Strain  Measure 

Recall  from  the  previous  sections  that  the  expression  for  the  nonlinear  strain  measure  is 


^5 


deg  =  In 


(4.36) 


Variational  Formulation 

The  functional  for  the  current  formulation  is  altered  such  that  the  bounds  of  the  integration 
are  performed  on  the  current  volume.  The  resulting  expression  for  the  internal  virtual  work 
now  becomes 


SYiint  =  dSs  as  Q.  I  =  51  Gg  (4-37) 

where  the  constraint  0,1  =  OqL  is  assumed  for  the  case  of  material  inelasticity^. 

^Recall  from  Chapter  2  the  volumetric  change  during  a  transformation  is  negligible. 
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Residual  and  Tangent  Arrays 


The  interpolation  functions  are  taken  to  be  the  same  as  the  linear  model,  resulting  in  the 
following  expression  for  the  variation  of  the  bar  length 


ndm  * 

(5/  =  ^  5Aui 

i=l 


I 


where  the  square  deformed  length  is  expressed  as 


(4.38) 


ndm  ndm 

1=1  i—l 


where  AXi  =  X2i  —  Xu  and  Aui  =  U2i  —  uu- 

Substitution  of  (4.38)  into  the  expression  for  the  internal  virtual  work  (4.37)  yields 


5n  =  5Aui 


ndfn  A 


i=l 


o  Q 


and  expressed  in  matrix  notation  we  have 


(5n  =  {  6uu 


(4.40) 


(4.41) 


where  Fi  =  (a  Q)^  is  the  internal  force. 

As  in  the  previous  section  we  introduce  one-way  thermal  effects  into  the  formulation  by 
adding  an  additional  nodal  quantity  or  degree  of  freedom,  namely  temperature  into  (4.41) 
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=  {  SUli  5Tu  Su2i  6T2i 

0 


-Fi  1 
0 
Fi 


>  . 


(4.42) 


The  vector  post-multiplying  the  virtual  displacements  and  temperatures  in  (4.42)  is  com¬ 
monly  referred  to  as  the  internal  load  vector  and  is  denoted  as  Pint  for  future  developments. 

To  facilitate  a  solution  of  the  mixed  boundary  value  problem  the  nonlinear  equation  (4.37) 
is  linearized  and  solved  by  Newton’s  method  as  a  sequence  of  linearized  problems.  Hence, 
linearizing  (4.37)  we  obtain 


d(OT)  =  d{5l)  as  ^  +  SI  das  ^  +  SI  as  dCl, 


(4.43) 


where 

d{Sl) 

dO, 
das 

and  Dm  and  Dt  are  the  mechanical  and  thermal  material  moduli  and  are  unchanged  from 
the  linear  case.  Substituting  (4.44)  into  (4.43)  with  the  aid  of  (4.38)  yields 

ndm  ndm 

d(m)  =  EE  SAui 

i=l  j=l 

and  expressed  in  matrix  notation  we  have 


+ 


Axi  Axj  {Dm  -  2a s) 


I  I 


I 


) 


I 


(4.45) 


ndm  ndm  -  ndm  a  ndm  * 

=  7  E  E  -  7  E  E  (4.44) 

i=l  j=l  2=1  j=l 

1 1  n^dyn  a 

”  ^  V  j2 

j=i 

=  Df^dSs  +  DxdT 
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d{5U)  = 


'  5uu  ' 

T 

kij  hi 

kij 

hi 

STu 

0  0 

0 

0 

^  5u2i 

> 

kij  hi 

kij 

hi 

ST2.  J 

0  0 

0 

0 

=  5t)^  k  dr] 


(4.46) 


where 


ki 


ij 


=  (^ 


^Xi  ^Xj  2c7^) 

~l  I  I 


Q  and 


^  Axj  ^  f2 


(4.47) 


The  matrix  post-multiplying  the  virtual  displacements  and  temperatures  in  (4.46)  is  com¬ 
monly  referred  to  as  the  tangent  matrix  and  is  denoted  as  fc.  As  presented  it  includes  both 
the  “material”  tangent  and  the  “geometric”  tangent. 


4.2.7  Solution  Procedure 

Noting  that  the  variations  6t)  in  the  previous  sections  are  arbitrary  we  obtain  the  finite 
element  residual  equations 


nelm 

Fintiv)  -  =  A  [F^ntiVe)  -  J^ext(^e)]  =  0  (4.48) 

e=i 


where  Pint  and  Fext  are  the  internal  and  external  load  vectors,  respectively,  e  is  an  element 
within  the  discretization,  nelm  denotes  the  total  number  of  elements  within  the  discretiza¬ 
tion,  r]e  =  {uii  Tji}  is  the  solution  vector  and  A  is  the  standard  finite  element  assembly 
operator. 

Linearizing  (4.48)  about  an  intermediate  state  yields 


L[Fint]  =  dve 


(4.49) 
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The  resulting  system  of  equations  involves  the  nodal  solution  vector  at  the  global  level 
expressed  as 


d-n  =  (4.50) 

where 


nelm 

=  A  (4.51) 

nelm 

=  A  \  pcxt  _  pint'}  (^) 

e~l  ^  Je  ' 

The  system  (4.50)  is  solved  and  then  the  unknown  fields  are  updated  additively  by 


,(Hi)  =  (4.52) 

The  process  is  repeated  within  a  particular  time  step  until  convergence  of  the  {k  +  1)** 
iterate  is  obtained,  the  solution  is  then  advanced  to  the  next  time  step  An  overview  of 
the  algorithmic  implementation  is  outlined  in  Table  4.1  below 

4.2.8  Numerical  Simulations 

We  consider  several  simulations  encompassing  both  isothermal  and  isostress  loading  histories 
to  demonstrate  the  features  of  the  model.  The  material  under  consideration  is  NiTi  and  its 
associated  material  parameters  can  be  found  in  Table  4.2. 

Isothermal  NiTi  Bar  Under  Cyclic  Loading,  T  <  T^/ 

We  consider  a  three  element  discretization  of  a  NiTi  bar  of  length  L  =  1  and  initial  area 
flo  =  1-  The  bar  is  fixed  against  translation  on  the  initial  end  {X  =  0)  and  loaded  at  the 
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_ Table  4.1:  Solution  outline  for  the  truss-bar  element 

Initialize  Variables 

Loop  over  time,  t 
Loop  until  convergence 
Initialize  Variables 
Loop  over  elements 
Compute  li  and  Axi 
Compute  axial  strain  Eg 
Compute  stress  and  moduli  ag,  Dm  and  Dt 
Compute  and  fc® 

Assemble  and 

End  loop 
Solve  K^'^^du  — 

End  loop 
End  loop 


Table  4.2:  Material  Properties  for  NiTi. (Brinson  k  Lammering  [1993]) 
Young’s  Moduli  Em  =  Ea  =  Q7  GPa 

Critical  stresses  for  de-twinning  =  100  MPa  and  =  170  MPa 
Martensite  production  temperatures  Tms  =  18.4  C  and  Tmf  =  9  C 
Austenite  production  temperatures  Tqs  =  34.5  C  and  Taf  =  49  C 
Austenite  production  slope  Cq  =  13.8  MPa/C 
Martensite  production  slope  Cm  =  8  MPa/C 
Maximum  transformation  strain  Ei  =  0.067 
Thermal  expansion  coefficient  a  =  6.5  /xstrain/C 


terminal  end  (X  =  L).  The  temperature  of  the  bar  is  held  constant  at  T  =  5  C,  which  is 
below  the  martensite  finish  temperature.  The  initial  state  is  assumed  to  be  100%  multiple 
variant  martensite  (i.e.  =  0.5)  and  the  loading  and  response  curves  are  outlined  in 

the  figures  below. 

From  Figures  4.10  and  4.11  we  can  conclude  that: 

1.  The  algorithm  accurately  determines  the  onset  and  completion  of  negative  and  positive 
martensitic  phase  transformations. 

2.  The  evolution  of  the  positive  and  negative  variants  is  continuous  and  smooth  during 
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Time  (sec) 


Tip  Displacement  (m) 


Figure  4.10:  Loading  and  response  curves  for  an  isothermal  NiTi  bar  under  cyclic  loading, 

T  <  Tmf. 

the  transformations.  Note,  the  production  of  one  variant  at  the  expense  of  the  other, 
while  their  sum  remains  unity,  i.e  100%  martensite. 

3.  Both  the  linear  and  non-linear  models  given  similar  results,  since  the  strain  level  is 
moderate. 
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Figure  4.11:  Time  history  of  the  internal  variables  for  an  isothermal  NiTi  bar  under  cyclic 
loading,  T  <Tmf- 

4.  For  the  completion  of  the  positive  variant  martensite  the  algorithm  has  a  possibility 
of  three  choices:  a)  elastic  unloading  with  b)  elastic  loading  =  1  and  c) 

positive  single  variant  martensitic  phase  transformation.  Upon  computing  the  three 
states  the  algorithm  determines  a)  is  inadmissible  since  da  >  0  and  =  0.  The 
remaining  paths  are  admissible  with  b)  selected  based  on  the  distance  measure  of  the 
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computed  and  previous  values  of  stress.  The  algorithm  continues  to  choose  b)  until 
the  current  time  step  converges.  Advancing  to  the  next  time  step  the  external  loading 
continues  to  increase  and  the  algorithm  determines  that  b)  is  the  only  admissible  path. 
A  similar  process  continues  throughout  the  remainder  of  the  loading  history. 


Isothermal  NiTi  Bar  Under  Cyclic  Loading,  Tas  <T  <Taj 

We  consider  a  three  element  discretization  of  a  NiTi  bar  of  length  L  =  1  and  initial  area 
f2o  =  1-  The  bar  is  fixed  against  translation  on  the  initial  end  (X  =  0)  and  loaded  at  the 
terminal  end  {X  =  L).  The  temperature  of  the  bar  is  held  constant  at  T  =  41.75  C,  which 
is  above  the  austenite  start  temperature  and  below  the  austenite  finish  temperature.  The 
initial  state  is  assumed  to  be  50%  multiple  variant  martensite  (i.e.  =  0.25)  and  the 

loading  and  response  curves  are  outlined  in  the  figures  below. 

From  Figures  4.12  and  4.13  we  can  conclude  that: 


1.  The  algorithm  accurately  determines  the  onset  and  completion  of  negative  and  positive 
martensitic  phase  transformations. 

2.  The  model  predicts  the  termination  of  a  austenite  transformation  when  crossing  over 
the  stress  axis. 

3.  The  evolution  of  the  positive  and  negative  variants  is  continuous  and  smooth  during 
the  transformations.  Note,  the  production  of  one  variant  at  the  expense  of  the  other. 

4.  Both  the  linear  and  non-linear  models  given  similar  results,  since  the  strain  level  is 
moderate. 

5.  For  the  termination  of  the  austenite  production  from  a  positive  stress  state  the  al¬ 
gorithm  has  a  possibility  of  three  choices:  a)  elastic  loading  with  ^  b)  elastic 

unloading  ^  =  ^apx  and  c)  an  austenitic  phase  transformation.  Upon  computing  the 
three  states  the  algorithm  determines  a)  is  inadmissible  since  dcr  <  0  and  ^  =  0.  The 
remaining  paths  are  admissible  with  c)  selected  based  on  the  distance  measure  of  the 
computed  and  previous  values  of  stress.  The  algorithm  continues  to  choose  c)  until 
the  current  time  step  converges.  Advancing  to  the  next  time  step  the  external  loading 
continues  to  decrease  and  the  only  admissible  path  determined  from  the  algorithm  is 
b). 
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Figure  4.12:  Loading  and  response  curves  for  an  isothermal  NiTi  bar  under  cyclic  loading, 

Tas<T  <  Taf. 

Isothermal  NiTi  Bzu*  Under  Cyclic  Loading,  T  >  To/ 

We  consider  a  three  element  discretization  of  a  NiTi  bar  of  length  L  =  1  and  initial  area 
Oo  =  1.  The  bar  is  fixed  against  translation  on  the  initial  end  {X  =  0)  and  loaded  at  the 
terminal  end  (X  =  L).  The  temperature  of  the  bar  is  held  constant  at  T  =  55  C,  which  is 
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Figure  4.13:  Time  history  of  the  internal  variables  for  an  isothermal  NiTi  bar  under  cyclic 
loading,  Tas  <  T  <  Taf. 

above  the  austenite  finish  temperature.  The  initial  state  is  assumed  to  be  100%  austenite 
(i.e.  =  0)  and  the  loading  and  response  curve  are  outlined  in  the  figures  below. 

Prom  Figures  4.14  and  4.15  we  can  conclude  that: 
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600| 


Figure  4.14:  Loading  and  response  curves  for  an  isothermal  NiTi  bar  under  cyclic  loading, 
T  >  Taf. 

1.  The  algorithm  accurately  models  the  pseudoelastic  behavior  by  determining  the  onset 
and  completion  of  negative  and  positive  martensitic  phase  transformations. 

2.  The  evolution  of  the  positive  and  negative  variants  is  continuous  and  smooth  during 
the  transformations.  Note,  only  one  variant  is  evolved  while  the  other  remains  zero, 
since  only  one  variant  can  exist  at  temperatures  greater  than  Taf. 
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Figure  4.15:  Time  history  of  the  internal  variables  for  an  isothermal  NiTi  bar  under  cyclic 
loading,  T  >Taf. 

3.  Both  the  linear  and  non-linear  models  given  similar  results,  since  the  strain  level  is 
moderate. 

4.  For  both  the  single  variant  martensite  and  austenite  productions  the  algorithmic  be¬ 
havior  is  similar  to  the  simulations  in  the  previous  sections. 
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5.  The  termination  of  the  single  variant  martensite  and  austenite  production  follow  from 
the  previous  cases. 


NiTi  Bar  Under  Thermo-mechanical  Cyclic  Loading 

We  consider  a  three  element  discretization  of  a  NiTi  bar  of  length  L  =  1  and  initial  area 
f2o  =  1.  The  bar  is  fixed  against  translation  on  the  initial  end  {X  =  0)  and  loaded  at  the 
terminal  end  {X  =  L).  The  temperature  of  the  bar  is  initially  held  at  65  C;  cooled  to  5  C 
under  zero  load;  then  loaded  and  unloaded  isothermally  at  5  C,  and  lastly  heated  to  65  C 
under  zero  load.  The  initial  state  is  assumed  to  be  100%  austenite  (i.e.  =  0)  and 

the  loading  and  response  curves  are  outlined  in  the  figures  below. 

From  Figures  4.16-4.18  we  can  conclude  that: 


1.  The  algorithm  accurately  models  the  shape  memory  effect  by  determining  the  onset 
and  completion  of  various  phase  transformations. 

2.  During  the  initial  cooling  process,  note  the  strain  developed  is  negligible  (thermal  ex¬ 
pansion  effects)  due  to  the  self-accommodating  state  which  arises  during  the  production 
of  multiple  variant  martensite. 

3.  During  the  loading  process  the  evolution  of  the  positive  variant  continues  to  grow  at 
the  expense  of  the  negative  variant  until  the  transformation  is  completed. 

4.  Upon  heating  the  bar  we  see  that  the  residual  strain  goes  to  zero  as  experimental 
observation  shows. 

5.  Both  the  linear  and  non-linear  models  given  similar  results,  since  the  strain  level  is 
moderate. 

6.  The  algorithmic  treatment  for  the  termination  of  the  single  variant  martensite  follows 
from  previous  sections.  Whereas,  for  the  austenite  transformation  the  algorithm  picks 
three  paths:  a)  elastic  loading  with  ^  =  ^„,  b)  an  austenitic  phase  transformation 
assuming  a  positive  value  for  the  stress  and  c)  an  austenitic  phase  transformation 
assuming  a  negative  value  for  the  stress.  Upon  computing  the  three  states  the  algorithm 
determines  a)  is  inadmissible  since  dT  >  0  and  ^  =  0.  The  remaining  paths  are 
admissible  with  c)  selected  based  on  the  distance  measure  of  the  computed  and  previous 
values  of  stress.  The  algorithm  continues  to  choose  c)  until  the  current  time  step 
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Figure  4.16:  Loading  and  response  curves  for  an  NiTi  bar  under  thermo-mechanical  cyclic 
loading 

converges.  The  transformation  advances  until  the  martensite  is  fully  depleted.  After 
the  martensite  fraction  is  depleted  the  algorithm  only  chooses  an  elastic  state. 
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Figure  4.17:  Time  history  of  the  internal  variables  for  an  NiTi  bar  under  thermo-mechanical 
cyclic  loading. 

NiTi  Truss-Bridge  Under  Cyclic  Loading 

We  consider  a  NiTi  cantilever  truss-bridge  consisting  of  53  elements  shown  in  Figure  4.19. 
The  bridge  is  20  units  long  by  2  units  deep  and  all  elements  have  an  initial  area  f2o  =  1- 
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Figure  4.18:  Thermo-mechanical  response  for  an  NiTi  bar  under  thermo-mechanical  cyclic 
loading. 


Figure  4.19:  Reference  configuration  for  the  truss-bridge. 

Four  simulations  are  performed  to  demonstrate  the  ability  of  the  algorithm  to  predict  the 
behavior  of  a  system  with  spatial  inhomogeneities  for  various  loading  conditions 

Simulation  #1 

The  truss-bridge  is  isothermally  loaded  to  a  peak  load  of  156  MN  and  then  unloaded  at 
temperature  of  5  C.  The  initial  state  of  the  material  is  100%  multiple  variant  martensite  (i.e. 
=  0.5).  The  response  curve  is  outlined  in  the  figure  below. 

Simulation  #2 

The  truss-bridge  is  isothermally  loaded  to  a  peak  load  of  25  MN  and  then  unloaded  at 
temperature  of  41.75  C.  The  initial  state  of  the  material  is  50%  multiple  variant  martensite 
(i.e.  =  ^0  =  0.25).  The  response  curve  is  outlined  in  the  figure  below. 
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Figure  4.20:  Response  curve  for  isothermal  conditions  at  5  C  for  a  NiTi  truss  bridge. 


Figure  4.21:  Response  curve  for  isothermal  conditions  at  41.75  C  for  a  NiTi  truss  bridge. 
Simulation  #3 

The  truss-bridge  is  isothermally  loaded  to  a  peak  load  of  35  MN  and  then  unloaded  at 
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temperature  of  55  C.  The  initial  state  of  the  material  is  100%  austenite  (i.e.  =  0). 

The  response  curve  is  outlined  in  the  figure  below. 


Tip  Displacement  (m) 


Figure  4.22:  Response  curve  for  isothermal  conditions  at  55  C  for  a  NiTi  truss  bridge. 
Simulation  #4 

The  truss-bridge  is  isothermally  loaded  to  a  peak  load  of  15  MN  and  then  unloaded  at 
temperature  of  5  C.  The  initial  state  of  the  material  is  100%  multiple  variant  martensite 
(i.e.  =  0.5).  The  truss-beam  is  then  heated  to  65  C  to  facilitate  the  recovery  of  the 

residual  strain  incurred  during  the  mechanical  loading  stage.  The  response  curve  is  outlined 
in  the  figure  below. 

From  Figures  4.20-4.23  we  can  conclude  that: 

1.  The  algorithm  accurately  models  the  global  effects  demonstrated  in  the  previous  uni¬ 
axial  cases. 

2.  For  temperatures  below  Taf  we  see  that  there  exists  a  global  evolution  in  the  response. 
This  evolution  is  manifested  as  a  result  of  the  fractions  changing  to  different  degrees 
from  the  initial  conditions  for  the  various  bars.  These  oscillations  reach  a  limit  cycle 
in  a  few  load  cycles  and  repeatability  of  the  response  is  then  seen. 
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Figure  4.23:  Thermo-mechanical  response  curve  for  a  NiTi  truss  bridge. 

3.  For  Figure  4.23  the  simulation  was  performed  using  the  linear  element  only  since  some 
bars  have  material  moduli  which  are  negative  (due  to  the  softening  effect  of  the  trans¬ 
formation)  resulting  in  an  ill  posed  problem  for  the  log-stretch  model.  An  alternative 
formulation  using  the  second  Piola-Kirchhoff  stress  and  the  Green-Lagrange  strain 
tensor  would  circumvent  this  problem. 

4.  The  difference  between  the  linear  and  non-linear  models  is  more  pronounced  since  the 
strains  are  larger  than  previous  cases  especially  for  elements  near  the  support. 


Summary  of  Numerical  Simulations 

We  have  considered  several  simulations  encompassing  both  isothermal,  isostress  and  non¬ 
isostress  paths  to  demonstrate  the  features  of  the  model  for  NiTi.  Several  noteworthy  char¬ 
acteristics  are: 


1.  The  algorithm  accurately  activates  the  correct  evolution  equation  for  the  corresponding 
state  on  the  phase  diagram. 

2.  The  algorithm  accurately  models  both  the  pseudoelastic  and  shape  memory  effect. 
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3.  The  algorithm  accurately  determines  the  onset  and  completion  of  the  various  phase 
transformations. 


4.  The  inclusion  of  a  simple  plastic  model  posses  no  difficultly  and  accurately  represents 
the  behavior  of  the  material,  see  Vandermeer  et.al.  [1981]  and  Govindjee  &  Kasper 
[1997]. 


4.3  Two-Dimensional  Beam  Element 


In  this  section  we  consider  the  formulation  of  a  two-dimensional  beam  element,  following  Simo 
et.  al.  [1984]  for  both  linear  and  nonlinear  kinematics;  see  Figure  4.24  for  a  representation  of 
the  geometry.  In  the  following  developments  we  regard  the  thermal  field  as  prescribed  such 
that  no  additional  degrees  of  freedom  are  necessary.  We  begin  with  a  brief  outline  of  the 
basic  notation  used  throughout  and  then  proceed  with  a  description  of  the  strain  measure 
used.  A  variational  equation  for  elastostatics  and  the  finite  element  interpolations  for  a  line 
element  are  reviewed.  Lastly,  the  finite  element  arrays  necessary  for  implementation  are 
developed. 


Figure  4.24:  Reference  and  deformed  configurations  for  the  beam  element. 


4.3.1  Notation 

We  consider  an  initially  straight  beam  of  length  L  and  cross  section  f2o  ^  with  smooth 
boundary  dQo  to  represent  a  bounded  reference  configuration  B  for  the  continuum  body.  We 
admit  the  decomposition  of  the  boundary  into  two  parts:  r.u;  C  dB  where  the  generalized 
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displacements  are  prescribed  as  w  =  w  and  C  dB  where  the  generalized  forces  are 
prescribed  as  R  =  R  subject  to 


dB  =  Ty,ll  Tr  and  F^  n  F/j  =  0 


(4.53) 


where  w  =  [u,  v,  6]^  denotes  the  generalized  displacements,  two  translations  and  a  rotations, 
and  R  =  [N,  V,  M]^  denotes  the  generalized  forces  of  the  beam,  a  normal  and  shear  force 
and  a  moment. 

For  subsequent  treatment  of  the  variational  formulation  we  distinguish  two  classes  of  func¬ 
tions  as  before,  namely,  the  space  of  admissible  solutions  and  the  space  of  variations. 

Let  U  be  the  space  of  admissible  generalized  displacements  written  as 


U  —  \  w  E  and  w  =  w  on  F^u}  (4.54) 

and  V  be  the  space  of  admissible  displacement  variations  written  as 

V  =  I  €  H\L)  and  5w  =  0on  F,^}  .  (4.55) 

4.3.2  Deformation  Measure 

We  assume  that  plane  sections  remain  plane,  but  not  necessarily  normal  after  deformation. 
Hence,  we  consider  the  two-dimensional  motion  of  the  centroid  of  the  beam  in  the  form 


a:i  =  +  m(Xi)  (4.56) 

X2  =  ■y(Xi)  . 


From  (4.56)  we  see  that  the  motion  of  the  line  of  centroids  of  the  beam  takes  the  form 
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Xo  =  +  ^2^2  —  (-^1  +  Cl  +  v{Xi)e2 


(4.57) 


By  application  of  (3.6)  the  deformation  gradient  of  the  line  of  centroids  for  the  motion  xo 
becomes 


FoEi  — 


dxo 

dx 


El  —  (1  +  +  v'e2 


(4.58) 


where  (•)'  denotes  differentiation  with  respect  to  the  Xi  coordinate  direction.  The  axial 
stretching  normal  to  the  deformed  beam  and  shearing  tangential  to  the  deformed  beam  are 
simply  the  dot  product  between  the  deformation  gradient  of  the  line  of  centroids  and  the 
bases  ga  on  the  deformed  configuration 


6i  =  FqEi  •  gi  =  (1  +  u')  cos{6)  +  v'sm{9)  (4.59) 

S2  =  FqEi  •  92  =  -  (1  +  u')  sin(^)  +  v'  cos{6)  . 

From  (4.59)  we  see  that  the  measures  of  deformation  take  the  form 


r  =  a¥- El 


r  <5i  -  n 

'  1  +  u'  ' 

62  >  where  F  =  < 

v' 

1  0'  J  1 

[  9'  J 

and  A  is  the  transformation  matrix  between  the  Gaussian  and  spatial  frames 


(4.60) 


Qi  =  Aci  with  A  = 


cos{6)  sin(0)  0 
—  sm{6)  cos(0)  0 
0  0  1 


Note  from  (4.60)  the  linearized  deformation  measures  becomes 


(4.61) 
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r  u' 

r  =  <  -e  +  v' 

[  9' 

which  agrees  with  Timoshenko  beam  theory,  see  Gere  &  Timoshenko  [1984]. 

Remark  4.3. 

If  the  beam  is  initially  rotated  within  the  reference  frame  an  appropriate  coordinate  trans¬ 
formation  is  performed  before  the  computation  of  the  deformation  measure,  see  Sack  [1989] 
for  further  details. 


4.3.3  Variational  Formulation 

An  approximate  solution  of  the  boundary  value  problem  is  constructed  from  a  variational 
statement  of  the  problem.  The  basic  field  equations  may  be  included  in  a  variational  state¬ 
ment  for  the  elasticity  problem  using  a  potential  functional.  Accordingly,  we  have  a  func¬ 
tional  in  which  the  deformation  vector,  w  eU  is  regarded  as  the  independent  variable.  The 
proposed  functional  11 :  ZY  — R  may  be  expressed  as 


n{w)  =  /  w{r)  dv  +  ne.t(tr)  (4.63) 

Jb 

where  W  (r)  is  a  stored  energy  function  and  for  conservative  external  loading 

Tlext{‘^)  =  ~  j  b  •  w  dV  —I  R  -  w  dA  (4.64) 

JB  Jtr 

b  being  the  body  force  per  unit  volume. 

We  may  state  the  problem  as:  Find  the  w  eU  which  makes  the  functional  n(tt7)  stationary 
for  all  admissible  variations  6w  G  V. 
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The  stationary  point  of  11  is  obtained  by  setting  to  zero  the  first  variation  of  (4.63)  with 
respect  to  the  independent  field.  Accordingly, 


<^n  =  J  ■  5r  dV  —  j  Sw  •  b  dV  —  j  5w  •  R  dA  =  0  in  B  (4.65) 

for  all  admissible  variations  Sw  €  V. 

We  may  recast  the  first  term  of  (4.65)  as 


(in,, 


int 


-L 


dW 

dr 


SrfddX 


(4.66) 


which  represents  the  internal  virtual  work,  where  for  a  general  singly  symmetric  cross  section 


dWjr) 

dr 


f  N 

0  —DmQ 

r  Ti  1 

0  Gfi  0 

< 

r2  i 

1  M  J 

~DmQ  0  DmI 

[raj 

=  Dr 


(4.67) 


where  fi,  Q  and  /  are  the  integrated  area,  first  moment  of  area  and  second  moment  of  area, 
respectfully  and  the  variation  of  the  deformation  is 


Sr  —  SA  F  +  (5F  =  Sd  A  F  +  A  (iF  . 


(4.68) 


Substituting  (4.68)  into  (4.65)  yields 


(in  =  y  (^se'^  iF^  +  (iie^  .1^)  ndx-  su^^t  ■ 


(4.69) 
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4.3.4  Finite  Element  Interpolations 


The  selection  of  the  generalized  displacement  interpolation  functions  for  the  two  node  initially 
straight  beam  element  will  be  the  same  as  those  for  the  bar  element  in  the  previous  sections, 
namely, 


s  s 

iVi  =  1  —  —  and  -/V2  =  y  . 

L/  Li 

The  position  vector  and  generalized  displacement  vector  follow  as  before. 


(4.70) 


4.3.5  Linear  Kinematic  Element 

We  proceed  by  defining  the  deformation  measures  used  for  the  formulation,  develop  the 
expression  for  the  internal  virtual  work  from  which  the  residual  equations  may  be  expressed. 
Lastly,  the  resulting  nonlinear  residual  equations  are  linearized  to  afford  a  solution  of  the 
nonlinear  set  of  equations  via  a  Newton’s  method. 


Deformation  Measure 

Recall  the  expression  for  the  linearized  deformation  measure  is 


u' 

-e^v' 

9' 


\ 


/ 


(4.71) 


where  u  is  the  displacement  in  the  Xi  coordinate  direction,  v  is  the  displacement  in  the  X2 
coordinate  direction  and  9  is  the  rotation  about  the  X3  coordinate  direction. 

We  may  express  the  virtual  generalized  deformation  vector,  with  the  aid  of  the  interpolation 
functions  and  (4.71)  as 
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6r  =  <  - 


r  6u'  Y 

5ui  ] 

T  p 

{  -so  +  ^ 

6vj 

1  J 

[  59j  J 

N'j 

0 

0 


0 

N'j 

0 


0 

0 


0 

0 


N'j  Nj 


-| 

■  1 

0 

0  ■ 

0 

1 

0 

0 

0 

1 

- 

0 

-1 

- 1 

o 

=  dvs^ 


(4.72) 


where  the  index  7  =  1, 2  is  the  local  node  number  and  vr  =  [uj^vj,  Oj]. 


Residual  and  Tangent  Arrays 

Substitution  of  (4.72)  into  the  expression  for  the  internal  virtual  work  yields 


SYiint  =  6w'^  j  b'^  rQ  dX  .  (4.73) 

The  vector  postmultiplying  the  virtual  generalized  displacements  in  (4.73)  is  commonly  re¬ 
ferred  to  as  the  internal  load  vector  and  is  denoted  as  Pint  for  future  developments. 

Remark  4.4. 

Recall  from  Chapter  3  that  the  stress  and  hence  the  resultant  stresses  R  are  in  general  a 
function  of  the  mechanical  and  thermal  state  and  internal  variables,  i.e.  R  =  R{e,T,^). 

To  facilitate  a  solution  of  the  mixed  boundary  value  problem  the  nonlinear  equation  (4.65) 
is  linearized  and  solved  by  a  Newton’s  method  as  a  sequence  of  linearized  problems.  Hence, 
linearizing  (4.65)  we  obtain 


^  d'^W{r)^ 


(4.74) 


where  for  convenience  define 
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d‘^W{r) 

5r2 


(4.75) 


and  substituting  (4.72)  and  (4.75)  into  (4.74)  and  collecting  terms  gives 

d{6n)  =  6w^  J  B^BB  dX  .  (4.76) 

The  matrix  postmultiplying  the  virtual  generalized  displacements  in  (4.76)  is  commonly 
referred  to  as  the  tangent  matrix  and  is  denoted  as  k. 


4.3.6  Non-Linear  Kinematic  Element 

We  proceed  by  defining  the  deformation  measures  used  for  the  formulation,  develop  the 
expression  for  the  internal  virtual  work  from  which  the  residual  equations  may  be  expressed. 
Lastly,  the  resulting  nonlinear  residual  equations  are  linearized  to  afford  a  solution  of  the 
nonlinear  set  of  equations  via  Newton’s  method. 


Deformation  Measure 

Recall  from  previous  sections  the  expression  for  the  nonlinear  strain  measure  is 


r  =  A¥-Ei  =  l  $2  >  .  (4.77) 

I  J 

where  (5i  —  1,  ^2  and  9'  are  interpreted  as  the  axial  strain,  shear  strain  and  curvature, 
respectively. 

We  may  express  the  virtual  generalized  deformation  vector  with  the  aid  of  the  interpolation 
functions  as 
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(4.78) 


'  6u'  " 
6v' 
69' 

69 

f  Sui  1 

-  N'j 

0 

0 

0 

6f]  =  < 

>  =1  Sv, 

0 

iv; 

0 

0 

1  J 

0 

0 

N'j 

Ni 

=  6w^ 


where  the  index  /  =  1, 2  is  the  local  node  number  and  vF  =  [uj,  vj,9i  . 


Residual  and  Tangent  Arrays 

Substitution  of  (4.78)  into  the  expression  for  the  internal  virtual  work  yields 


(illrt,  =  Sv>‘'  J  RQ  dX 


(4.79) 


where 


■  ■ 

—  sin(0) 

cos(0)  0 

and  A  = 

—  cos(0) 

0 

—  sin(0)  0 

0  0 

(4.80) 


The  vector  postmultiplying  the  virtual  generalized  displacements  in  (4.79)  is  commonly  re¬ 
ferred  to  as  the  internal  load  vector  and  is  denoted  as  Pint  for  future  developments. 

To  facilitate  a  solution  of  the  mixed  boundary  value  problem  the  nonlinear  equation  (4.65) 
is  linearized  and  solved  by  Newton’s  method  as  a  sequence  of  linearized  problems.  Hence, 
linearizing  (4.65)  we  obtain 


d^Wjr) 


Sr  + 


dW{r) 

dr 


d{6r)  dX 


(4.81) 


where 
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(i{^6r'j  —  +  (JilcnF  +  ciA(5]F 


(4.82) 


d{5A)  =  69Ad9  and  A  — 


-  cos{9)  -  sin(0)  0 
sin(0)  —  cos{9)  0 
0  0  0 


(4.83) 


Substituting  (4.82)  and  (4.83)  into  (4.81)  and  collecting  terms  gives 


d{5n)  =  5w^ 


dX  + 


J  B^Gb  dX^ 


dw 


(4.84) 


where 


0  A^R 
if  A  R 


(4.85) 


The  matrix  postmultiplying  the  virtual  generalized  displacements  in  (4.84)  is  commonly 
referred  to  as  the  tangent  matrix  and  is  denoted  as  fc.  Note  the  first  term  of  (4.84)  represents 
the  material  part  of  the  tangent  and  the  second  term  the  geometric  part  of  the  tangent. 


4.3.7  Solution  Procedure 

Noting  that  the  variations  dr}  in  the  previous  sections  are  arbitrary  we  obtain  the  finite 
element  residual  equations 


nelm 

Fintil)  -  FeM  =  A  [Fl^tiVe)  "  ■Fext(^e)]  =  0  (4.86) 

e=l 
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where  Fint  and  Fext  are  the  internal  and  external  load  vectors,  respectively,  e  is  an  element 
within  the  discretization,  nelm  denotes  the  total  number  of  elements  within  the  discretiza¬ 
tion,  T}e  =  {uii  Vii  9ii}  is  the  solution  vector  and  A  is  the  standard  finite  element  assembly 
operator. 

Linearizing  (4.86)  about  an  intermediate  state  yields 


L[Fint]  =  dve  (4.87) 

The  resulting  system  of  equations  involves  the  nodal  solution  vector  at  the  global  level 
expressed  as 

dr)  =  (4.88) 

where 


nelm 

=  A  [fc]f )  (4.89) 

e=l  ^ 
nelm 

=  A  (  pext  _ 

6=1  ^  -le  ' 

The  system  (4.88)  is  solved  and  then  the  unknown  fields  are  updated  additively  by 


(4.90) 

The  process  is  repeated  within  a  particular  time  step  until  convergence  of  the  {k  +  1)*^ 
iterate  is  obtained,  the  solution  is  then  advanced  to  the  next  time  step  t„+i.  An  overview  of 
the  algorithmic  implementation  is  outlined  in  Table  4.3  below 
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_ Table  4.3:  Solution  outline  for  the  beam  element 

Initialize  Variables 

Loop  over  time,  t 

Loop  until  convergence 
Initialize  Variables 
Loop  over  elements 
Compute  Ni  and 

Loop  over  integration  points  along  length 
Compute  r 

Loop  over  integration  points  thru  depth 
Compute  axial  strain  e* 

Compute  stress  and  material  moduli  ag  and  Dm 
Compute  Q,  Q  and  I 
Compute  and  accumulate  R 
End  loop 

Compute  Pint  and  k 
End  loop 

Assemble  and 
End  loop 

Solve  =  R^'^^ 

End  loop 
End  loop 


4.3.8  Numerical  Simulations 

We  consider  several  simulations  encompassing  isothermal  and  non-isostress  loading  histories 
to  demonstrate  the  features  of  the  model.  The  material  under  consideration  is  NiTi  and  its 
associated  material  parameters  can  be  found  in  Table  4.2 


Isothermal  NiTi  Beam  Under  Cyclic  Loading,  T<Tmf 

We  consider  a  ten  element  discretization  of  a  NiTi  cantilever  beam  of  length  L  —  20m  having 
a  rectangular  cross  section  of  2m  x  Im.  The  cross  section  is  divided  into  four  layers,  each 
of  which  is  evaluated  using  a  5-pt  Gauss-Lobatto  quadrature  rule.  Standard  1-pt  Gaussian 
quadrature  was  used  along  the  axis  of  each  element.  The  beam  is  fixed  against  translations 
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and  rotations  on  the  initial  end  (X  =  0)  and  loaded  via  displacement  control  in  the  X2 
direction  at  the  terminal  end  {X  =  L).  The  temperature  of  the  bar  is  held  constant  at 
T  =  5  C,  which  is  below  the  martensite  finish  temperature.  The  initial  state  is  assumed 
to  be  100%  multiple  variant  martensite  (i.e.  =  0.5)  and  the  loading  and  response 

curves  are  outlined  in  the  figures  below. 


Tip  Displacement  (m) 


Figure  4.25:  Loading  and  response  curves  for  an  NiTi  beam  under  a  under  cyclic  loading, 
T  <  Tmf- 
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From  Figure  4.25  we  can  make  the  following  observations: 


1.  The  algorithm  accurately  determines  the  onset  and  completion  of  negative  and  positive 
martensitic  phase  transformations  for  each  quadrature  point  through  the  depth  of  the 
beam  resulting  in  a  behavior  which  is  qualitatively  correct. 

2.  Although  the  strain  levels  within  the  beam  are  moderate  the  rotations  are  large  and 
their  effects  are  pronounced. 


Isothermal  NiTi  Beam  Under  Cyclic  Loading,  Tas  <T  <Taf 

We  consider  a  ten  element  discretization  of  a  NiTi  cantilever  beam  of  length  L  =  20m 
having  a  rectangular  cross  section  of  2m  x  Im.  The  cross  section  is  divided  into  four  layers, 
each  of  which  is  evaluated  using  a  5-pt  Gauss-Lobatto  quadrature  rule.  Standard  1-pt 
Gaussian  quadrature  was  used  along  the  axis  of  each  element.  The  beam  is  fixed  against 
translations  and  rotations  on  the  initial  end  (A  =  0)  and  loaded  via  displacement  control  in 
the  X2  direction  at  the  terminal  end  {X  =  L).  The  temperature  of  the  bar  is  held  constant 
at  T  =  41.75  C,  which  is  above  the  austenite  start  temperature  and  below  the  austenite 
finish  temperature.  The  initial  state  is  assumed  to  be  50%  multiple  variant  martensite  (i.e. 

=  0.25)  and  the  loading  and  response  curves  are  outlined  in  the  figures  below. 

From  Figure  4.26  we  can  make  the  following  observations: 

1.  The  algorithm  accurately  determines  the  onset  and  completion  of  negative  and  positive 
martensitic  phase  transformations  for  each  quadrature  point  through  the  depth  of  the 
beam  resulting  in  a  behavior  which  is  qualitatively  correct. 

2.  The  model  predicts  the  termination  of  a  austenite  transformation  when  the  net  reaction 
is  zero. 

3.  Although  the  strain  levels  within  the  beam  are  moderate  the  rotations  are  large  and 
their  effects  are  pronounced. 


Isothermal  NiTi  Beam  Under  Cyclic  Loading,  T  >  Taf 

We  consider  a  ten  element  discretization  of  a  NiTi  cantilever  beam  of  length  L  —  20m  having 
a  rectangular  cross  section  of  2m  x  Im.  The  cross  section  is  divided  into  four  layers,  each 
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Figure  4.26:  Loading  and  response  curves  for  an  NiTi  beam  under  a  under  cyclic  loading, 
Tas<T  <  Taf. 


of  which  is  evaluated  using  a  5-pt  Gauss-Lobatto  quadrature  rule.  Standard  1-pt  Gaussian 
quadrature  was  used  along  the  axis  of  each  element.  The  beam  is  fixed  against  translations 
and  rotations  on  the  initial  end  (X  =  0)  and  loaded  via  displacement  control  in  the  X2 
direction  at  the  terminal  end  {X  =  L).  The  temperature  of  the  bar  is  held  constant  at 
T  =  55  C,  which  is  above  the  austenite  finish  temperature.  The  initial  state  is  assumed  to 
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be  100%  austenite  (i.e.  =  0)  and  the  loading  and  response  curve  are  outlined  in  the 

figures  below. 


Figure  4.27:  Loading  and  response  curves  for  an  NiTi  beam  under  a  under  cyclic  loading, 
T  >  T,f. 

From  Figure  4.27  we  can  make  the  following  observations: 

1.  The  algorithm  accurately  determines  the  onset  and  completion  of  negative  and  positive 
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martensitic  phase  transformations  for  each  quadrature  point  through  the  depth  of  the 
beam  resulting  in  a  behavior  which  is  qualitatively  correct. 

2.  From  the  experimental  data  within  Auricchio  et.al.  [1995]  the  model  qualitatively 
captures  the  essential  behavior  of  the  pseudoelastic  effect. 

3.  Although  the  strain  levels  within  the  beam  are  moderate  the  rotations  are  large  and 
their  effects  are  pronounced. 


NiTi  Beam  Under  Thermo-mechanical  Cyclic  Loading 

We  consider  a  ten  element  discretization  of  a  NiTi  cantilever  beam  of  length  L  =  20m  having 
a  rectangular  cross  section  of  2m  x  Im.  The  cross  section  is  divided  into  two  layers,  each 
of  which  is  evaluated  using  a  3-pt  Gauss-Lobatto  quadrature  rule.  Standard  1-pt  Gaussian 
quadrature  was  used  along  the  axis  of  each  element.  The  beam  is  fixed  against  translations 
and  rotations  on  the  initial  end  {X  =  0)  and  loaded  via  load  control  in  the  X2  direction  at 
the  terminal  end  (A  =  L).  The  bar  is  initially  loaded  and  unloaded  isothermally  at  5  C, 
and  lastly  heated  to  80  C  under  zero  load.  The  initial  state  is  assumed  to  be  100%  multiple 
variant  martensite  (i.e.  ^0  =  ^0  —  0-5)  and  the  loading  and  response  curves  are  outlined  in 
the  figures  below. 


Figure  4.28:  Loading  and  response  curves  for  an  NiTi  beam  under  a  under  cyclic  loading. 


2-100 


Figure  4.29:  Thermo-mechanical  response  curves  for  an  NiTi  beam  under  a  under  cyclic 
loading. 

From  Figures  4.28  and  4.29  we  can  conclude  that: 


1.  The  algorithm  accurately  models  the  shape  memory  effect  by  determining  the  onset 
and  completion  of  various  phase  transformations  for  each  quadrature  point. 

2.  Upon  heating  the  bar  we  see  that  the  residual  deformation  goes  to  zero  as  experimental 
observation  shows. 

Remairk  4.5. 

The  simulations  were  also  performed  using  various  quadrature  orders  and  total  number  of 
layers.  Quadrature  rules  greater  than  3  resulted  in  unconverged  states  during  the  heating 
cycle,  specifically  during  the  production  of  austenite.  Initial  findings  indicate  that  since 
the  resultant  force  on  the  cross-section  is  zero  multiple  stress  distribution  are  admissible, 
thereby  making  state  determination  indeterminate.  Current  work  in  underway  to  determine 
alternative  schemes  for  the  determination  of  the  state  which  maximizes  the  dissipative  energy. 
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4.3.9  Summary  of  Numerical  Simulations 

We  have  considered  several  simulations  encompassing  isothermal  load  histories  to  demon¬ 
strate  the  features  of  the  model  for  NiTi,  in  the  setting  of  a  linear  and  non-linear  two 
dimensional  beam.  Several  noteworthy  characteristics  are: 


1.  The  algorithm  accurately  activates  the  correct  evolution  equation  for  the  corresponding 
state  on  the  phase  diagram. 

2.  The  algorithm  accurately  models  both  the  pseudoelastic  and  shape  memory  effect. 

3.  The  algorithm  accurately  determines  the  onset  and  completion  of  the  various  phase 
transformations. 

4.  Since  the  beam  has  multiple  quadrature  through  the  depth  it  provides  a  rigorous 
examination  of  the  algorithm,  due  to  the  possibility  of  complex  stress  states. 


4.4  Summary 


In  the  previous  sections  we  have  consider  the  formulation  of  a  thermomechanical  constitutive 
model  for  shape  memory  alloys  and  its  implementation  into  a  multi-dimensional  truss-bar  and 
two-dimensional  beam  element.  The  simulations  presented  in  this  work  provide  a  valuable 
tool  for  the  analysis  and  design  of  shape  memory  devices.  Some  of  the  key  features  emanating 
from  this  work  are  outlined  below: 

1.  The  development  of  a  new  constitutive  model  which  accounts  for  both  compressive  and 
tensile  states  of  stress  and  the  associated  variant  productions  is  presented. 

2.  We  introduce  a  new  algorithm  for  the  determination  of  the  state  which  accounts  for 
nested  elastic  and  inelastic  regimes  based  upon  the  use  of  modified  trial  state  variables. 

3.  The  ability  of  the  model  to  replicate  the  quantitative  behavior  for  the  shape  memory 
and  pseudoelastic  effect  is  presented  for  both  truss  and  beam  elements  for  linear  and 
finite  kinematics. 

4.  The  algorithm  developed  for  state  determination  accurately  computes  onset,  evolution 
and  completion  of  the  various  transformations  which  may  occur. 
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Part  III 


A  Nonlinear  Composite  Shell  Theory 


1.  GEOMETRICALLY  NONLINEAR  APPROACHES 


FOR  FINITE  ELEMENT  ANALYSIS 

1.1  Introduction 

A  general  purpose  finite  element  code  having  powerful  and  reliable  geometric  and  mate¬ 
rial  nonlinear  capabilities  is  an  essential  tool  for  today’s  structural  engineer.  In  this  work  we 
have  endeavored  to  extend  a  multi-layered  shear  deformable  composite  shell  theory,  initially 
designed  for  linear  analysis  only,  to  the  nonlinear  regime.  In  the  first  phase  of  the  enhance¬ 
ment  for  nonlinear  analysis,  a  geometrically  nonlinear  capability  for  large  displacements  and 
rotations  but  small  strains  has  been  implemented. 

A  comprehensive  literature  review  was  conducted  on  geometrically  nonlinear  finite  ele¬ 
ment  approaches  and  nonlinear  solution  techniques  developed  in  the  last  two  decades.  As 
discussed  in  the  following  sections,  the  element  co-rotational  procedure  was  selected  along 
with  a  nonlinear  solution  technique  for  handling  the  geometric  nonlinearity.  This  co-rotational 
procedure  is  ideally  suited  to  extending  an  existing  linear  finite  element  code  to  perform  ge¬ 
ometrically  nonlinear  analysis.  The  theoretical  issues  related  to  the  implementation  of  the 
co-rotational  procedure  are  discussed  and  the  close  relationship  between  the  updated  La- 
grangian  formulation  and  the  co-rotational  procedure  is  revealed.  All  the  approximations 
made  in  the  co-rotational  procedure  as  well  as  the  restrictions  and  limitations  caused  by 
them  are  also  identified. 

1.2  Background  on  Total  and  Updated  Lagrangian  Formulation 

In  the  past  two  decades,  significant  advancements  have  been  made  in  the  development  of 
efficient  and  convenient  finite  element  theories  for  large  deflection  and  inelastic  response  of 
complex  structures.  Several  formulation  strategies  and  procedures  axe  available  for  accommo¬ 
dating  large  rotation  capabilities.  The  advances  made  have  been  prompted  by  an  improved 
understanding  of  the  physical  principles  involved  and  the  continuing  advances  in  computer 
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technology. 


When  large  displacements,  large  rotations  and  even  finite  strain  nonlinear  constitutive 
relations  are  involved  in  problems,  an  incremental  solution  technique  is  generally  required  and 
either  a  total  Lagrangian  or  an  Eulerain  coordinate  system  is  employed  to  establish  the  finite 
element  equilibrium  equations.  The  formar  coordinate  system  results  in  equilibrium  equations 
in  terms  of  the  initial  undeformed  configuration  whereas  the  latter  coordinate  system  results  in 
equilibrium  equation  with  respect  to  the  final  deformed  configuration.  However,  because  the 
final  deformed  configuration  is  determined  as  a  part  of  the  solution,  an  Eulerian  formulation 
reduces  to  an  updated  Lagrangian  formulation,  where  the  equifibrium  state  obtained  for 
the  previous  load  increment  (or  iteration)  is  employed  as  the  reference  configuration  for  the 
current  load  increment  (or  iteration). 


1.3  Basic  Philosophy  of  the  Co-rotational  Procedure 

Unfortunately,  neither  the  total  nor  the  updated  Lagrangian  formulation  can  be  simply 
incorporated  into  existing  linear  finite  element  codes.  It  is  well  known  in  continuum  mechanics 
that  the  motion  of  a  continuous  medium  can  always  be  decomposed  into  a  rigid  body  motion 
followed  by  a  pure  deformation.  If  the  finite  element  is  small  enough  to  provide  a  valid 
approximation  of  the  continuous  body,  then  the  motion  of  each  individual  finite  element  can 
be  similarly  decomposed.  In  particular,  the  motion  of  a  flexible  shell  or  beam  finite  element 
can  be  conceived  as  a  finite  rigid  body  motion  followed  by  relatively  small  pure  deformation. 
The  latter  is  not  involved  with  geometrical  nonlinearities  and  the  former  is  not  involved  in  the 
constitutive  equations  for  the  element.  If  the  rigid  body  motion  part  is  efiminated  from  the 
total  displacements,  the  deformational  part  of  the  motion  is  always  a  small  quantity  relative 
to  the  local  element  axes.  Since  the  strain  is  small,  the  deformational  displacement  derivatives 
with  respect  to  the  local  element  coordinates  (convected  coordinates)  must  be  in  the  same 
order  of  magnitude  as  the  strains.  As  a  result,  the  second  order  terms  in  the  Green-Lagrange 
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strain  formulation  can  be  dropped  [3]  and  the  linear  finite  element  theory  can  be  adopted  for 
geometrically  nonlinear  analysis. 

Based  on  this  concept,  Belytschko  and  Hsieh  [4]  developed  a  co-rotational  procedure  for 
nonlinear  transient  analysis  and  Wempner  [5]  proposed  a  finite  rotation  and  small  strain  finite 
element  theory  for  flexible  shell  elements  in  static  analysis.  In  [4],  each  element  is  associated 
with  a  local  convected  coordinate  system  which  is  assumed  to  rotate  relative  to  the  fixed 
global  coordinate  system  by  the  rigid  body  rotation  of  that  element.  With  this  convected 
coordinate  system,  the  total  displacements  and  rotations  were  easily  decomposed  and  the 
geometrical  nonlinearities  caused  by  large  rotations  were  treated  entirely  by  transformations 
between  the  global  and  the  convected  coordinates. 

Since  the  co-rotaitonal  procedure  reported  in  [4,5],  many  elements  based  on  this  proce¬ 
dure  have  been  developed  for  geometrically  nonlinear  analysis.  Harrigmoe  and  Bergan  [6] 
successfully  applied  the  co-rotational  procedure  to  linear  fiat  triangular  and  rectangular  shell 
elements  for  large  rotation  analysis.  In  their  study,  the  local  convected  coordinate  system 
was  natmally  defined  by  using  the  direction  cosines  of  the  mid-surface  normal  and  one  edge 
of  the  element.  Excellent  results  for  several  typical  problems  involving  large  rotations  and 
buckling  and/or  bifurcation  points  were  obtained. 


1.4  Solution  Techniques  for  Nonlinear  Finite  Element  Equations 

The  nonlinear  finite  element  equations  of  motion  are  normally  expressed  as  a  set  of  non¬ 
linear  algebraic  equations  containing  nodal  displacements  and  external  load  vectors.  Among 
them,  each  algebraic  equation  states  the  equifibrium  condition  of  an  individual  degree-of- 
ffeedom  in  the  finite  element  model.  For  the  case  of  proportional  loading,  the  external  load 
vector  can  always  be  expressed  as  the  product  of  an  arbitrarily  defined  reference  external 
load  vector  and  a  scalar  called  the  load  parameter.  Actually  the  load-displacement  curves 
are  determined  by  these  N  nonlinear  algebraic  equations  about  {N  -|- 1)  unknowns  where  N 
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is  the  total  degrees-of-freedom  of  the  finite  element  model. 

The  finite  element  equations  can  be  solved  in  various  ways,  for  instance,  by  the  pure 
incremental  method  (without  equilibrium  iteration),  the  iterative  method  based  on 
the  secant  stiffness  relations  or  the  combined  incremental-iterative  method  based 
on  the  Newton-Raphson  scheme.  The  finite  element  stabihty  analysis  is  characterized  by  a 
singular  or  ill-conditioned  tangent  stiffness  matrix  at  or  in  the  vicinity  of  the  critical  load 
point.  In  the  full  Newton-Raphson  method  where  the  stiffness  matrix  is  modified  at  the 
end  of  each  iteration,  it  is  quite  possible  that  the  critical  point  may  be  hit  exactly  and  as  a 
result  the  iterative  procedure  break  down.  However,  the  modified  Newton-Raphson  method  is 
better  at  avoiding  the  critical  point  and  is  therefore  highly  recommended  in  buckling  analysis 
[13].  According  to  several  researchers  [12-15],  and  our  experiences,  the  exact  singularity  in 
tangent  stiffness  matrix  is  extremely  unlikely  to  be  encountered  in  practical  computation  if  the 
modified  Newton-Raphson  method  is  employed.  Therefore  we  employ  the  Modified  Newton 
Raphson  Method  together  with  the  hne-search  technique  to  solve  the  resulting  nonlinear 
system  of  equations. 
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2.  ELEMENT  COROTATIONAL  FORMULATION 


FOR  GEOMETRICALLY  NONLINEAR  ANALYSIS 
OF  MULTILAYERED  SHELLS 

In  this  section  the  theoretical  formulation  of  the  element  co-rotational  procedure  for  geo¬ 
metrically  nonhnear  analysis  is  presented.  Following  a  brief  introduction  identifying  the  main 
attributes  of  the  total  and  updated  Lagrangian  formulations,  the  co-rotational  procedure  is 
derived  from  the  updated  Lagrangian  method  consistently  by  making  a  series  of  approxima¬ 
tions  based  on  the  assumption  of  small  strain.  The  sources  of  possible  numerical  errors  and 
restrictions  on  the  co-rotational  procedure  resulting  from  the  approximations  in  its  derivation 
are  also  discussed. 

To  present  the  nonlinear  finite  element  theory,  we  consider  the  motion  of  a  general  body  in 
a  stationary  Cartesian  coordinate  system.  It  is  assmned  that  the  general  body  can  experience 
large  displacements  and  rotations,  large  strains,  and  a  nonlinear  constitutive  response.  In 
static  analysis,  the  time  variable  is  employed  just  for  convenience  in  describing  the  loading 
process.  It  is  obvious  that  the  objective  of  nonlinear  analysis  is  to  evaluate  the  equilibrium 
positions  of  the  complete  body  at  a  series  of  discrete  time  points  0,  At,  2At,  3At...,  where  the 
time  increment  At  corresponds  to  the  load  increment.  In  an  incremental  solution  procedure, 
assuming  that  the  solutions  for  the  static  and  kinematic  variables  for  all  the  time  points  from 
0  to  t  have  been  obtained,  we  have  to  develop  a  solution  strategy  to  solve  for  the  next  required 
equilibrium  position  corresponding  to  time  t  H-  At. 

In  the  Lagrangian  incremental  analysis  approach,  we  can  express  the  equihbrium  of  the 
body  at  time  t  +  At  using  the  principle  of  virtual  work.  Using  tensor  notation,  this  principle 
requires  that 

/  <Tn+l  <5en+i  dfl  =  (1) 

J 

where  CTn+i  is  the  Cauchy  stress  tensor,  is  the  infinitesimal  strain  tensor,  and  the  6 
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represents  variation.  Both  <t„+i  and  are  referred  to  the  configuration  at  time  t  +  At. 
The  external  virtual  work,  fn^i  can  be  expressed  as 

/;?!=  /  (2) 

jQt  JTt 

where  fn+i  a^nd  fn+i  indicate  the  externally  applied  body  and  surface  force  vectors,  respec¬ 
tively,  and  6u  is  the  virtual  displacement  vector.  Q,t  and  F*  denote  the  volume  and  surface 
corresponding  to  the  equilibrium  configuration  at  time  t  +  At. 

It  should  be  noted  that  the  principle  of  virtual  work  presented  in  equation  (1)  is  established 
in  the  unknown  equilibrium  configuration  at  time  t  +  At.  This  is  a  fundamental  difference 
compared  with  linear  analysis  in  which  it  is  assumed  that  the  displacements  are  infinitesimally 
small  so  that  the  configuration  of  the  body  does  not  change.  To  overcome  this  difficulty,  either 
the  total  Lagrangian  or  the  updated  Lagrangian  formulation  is  usually  employed. 

Remark:  In  the  following  we  will  follow  the  notation  from  Marsden  and  Hughes  [1983]. 


2.1  Total  Lagrangian  Formulation 

In  the  total  Lagrangian  formulation,  the  principle  of  virtual  work  (1)  is  estabhshed  with 
respect  to  the  initial  undeformed  configuration  at  time  0. 

[  Sn+i  5En+l  d^o  =  F^l\  (3) 

Jno 

where  En+i  is  the  Green-Lagrange  strain  tensor  given  by 


En+l  =  UK,j)n+l 


(4) 


In  (4),  the  displacement  gradients  Ui,j  with  respect  to  the  initial  undeformed  configuration 
at  time  0  are  defined  as 


dUj 

dXj 


(5) 


where  U  indicates  the  total  displacements  referred  to  the  initial  configuration  and  X  repre¬ 
sents  the  coordinates  in  the  initial  configuration.  Sn+i  denotes  the  second  Piola-Kirchhoff 
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stress  tensor  which  is  related  to  the  Cauchy  stress  tensor  by  a  stress  transformation  as  follows: 


^n+l  —  ^n+l  -^n+1 


(6) 


where  J  is  the  Jacobian  of  transformation  and  F  ^  is  the  inverse  of  the  deformation  gradient 
matrix,  namely, 


_  dxi  1  _  ^ 

~dXi  ’  "  dxi 


(7) 


Because  the  second  Piola-Kirchhoff  stress  tensor  and  the  Green-Lagrange  strain  tensor  are 
energetically  conjugate,  they  can  be  used  as  a  pair  in  the  principle  of  virtual  work. 


Remark:  It  should  be  noted  that  the  Cauchy  stress  tensor,  which  corresponds  to  the  internal 
force  per  unit  area  in  the  final  deformed  configuration,  represents  the  true  stress  state  in  the 
deformed  body.  However,  the  second  Piola-Kirchhoff  stress  tensor,  which  is  interpreted  as 
internal  force  per  unit  area  in  the  initial  configuration,  does  not  represent  the  true  stress  state 
in  the  structure  and  cannot  be  used  directly  in  engineering  design. 


Defining  the  displacement  increments  within  the  load  step  from  t  to  f  -h  At  as  follows: 


AU  =  Un+l  -  Un 


(8) 


We  can  divide  the  Green-Lagrange  strain  increment  AE  =  ^/n+i  ~  Fn  into  linear  and  non¬ 
linear  parts  as 

AE  =  E^  +  E^’'  (9) 


where 

E\j  =  ^{Uij  +  Ujj  +  {UK,l)nUK,J  +  UkJ  {UK,j)n)  (10) 

and 

=  \Uk.,Ukj  (11) 

Remark:  The  initial  displacement  terms  {UK,i)n  and  {UK,j)n  introduce  the  initial  displace¬ 
ment  effect  in  the  elastic  stiffness  matrix  and  the  internal  force  vector. 
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The  second  Piola-Kirchhoff  stress  tensor  in  equation  (3)  can  be  expressed  as  the  summa¬ 
tion  of  the  stress  components  at  time  t  and  a  small  increment  as 

Sn+l  =  Sn  +  AS  (12) 

Substituting  (9)  and  (12)  the  principle  of  virtual  work  (3),  denoting  6En+i  =  6AE  and 
neglecting  the  quadratic  terms,  we  obtain  the  linearized  equations  of  motion  in  incremental 
form 

[  E^:C:6E^dn+  f  SndE'^^dn  =  F^l\-  f  SnSE^dQ  (13) 

Juo  JUo 

where  the  constitutive  tensor  C  is  used  to  relate  the  incremental  second  Piola-Kirchhoff  stress 
to  the  linear  part  of  the  incremental  Green-Lagrange  strain  in  the  first  term  on  the  left  hand 
side  in  (13)  as 

AS  =  C:E^  (14) 

Remark:  It  should  be  pointed  out  that  because  constitutive  equation  (14)  is  given  in 

incremental  form,  is  can  be  understood  as  a  linearization  of  any  general  nonlinear  constitutive 
relation. 

The  linearized  principle  of  virtual  work  presented  in  (13)  can  be  employed  to  calculate 
an  increment  in  the  displacements  for  prescribed  external  load  at  t  At,  which  is  then  used 
to  evaluate  the  approximation  to  the  total  displacements,  strains  and  stresses  corresponding 
to  the  time  t  +  At.  The  displacement  approximations  corresponding  to  t  +  At  are  simply 
obtained  by  adding  the  calculated  displacement  increments  to  the  displacements  at  time  t. 
The  strain  and  stress  approximations  are  then  evaluated  from  the  displacements  using  the 
kinematic  relation  and  the  constitutive  relation. 

After  obtaining  the  approximate  displacements,  strains  and  stresses,  we  check  the  differ¬ 
ence  between  external  virtual  work  and  the  internal  virtual  work  evaluated  from  the  obtained 
approximate  static  and  kinematic  variables  at  time  t  +  At.  Denoting  the  approximate  val¬ 
ues  with  a  superscript  {i)  in  anticipation  that  an  iteration  will,  in  general,  be  necessary,  the 
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residual  due  to  linearization  (13)  is 

Residiud  =  -  /  Sj,‘l,SE^l,dQ  (15) 

J  f2o 

Remark:  It  is  important  to  note  that  the  right  hand  side  of  (15)  is  equivalent  of  the  right 
hand  side  of  (13).  However,  in  (15)  the  current  value  of  stress  and  strain  variables  are 
employed.  The  above  considerations  show  that  the  right  hand  side  terms  in  (13)  and  (15)  are 
the  “out-of-balance  virtual  work” . 


In  order  to  reduce  the  out-of-balance  virtual  work  to  within  a  certain  convergence  measure, 
iterations  are  normally  required.  Denoting  the  iteration  number  by  superscript  (i),  we  have, 
for  i=  1,2, 3, ... 


E‘^*”  :  C:  dn  +  /  S„+i 


se:‘T  =  -f;?!  -  / 


n+1 


si+l  SEihda 


(16) 


Dming  the  iterations  the  displacements  are  updated  as  follows 


n+1 


"n+1 


(17) 


where  in  (17)  indicates  the  displacement  increment  obtained  from  iteration  i  to  i-|-  1. 

The  initial  displacement  in  any  typical  step  is  the  last  converged  displacement  in  the  previous 
step  and  is  given  by 

C/it  =  U„  (18) 


2.2  Updated  Lagrangian  Formulation 

The  updated  Lagrangian  formulation  also  can  be  derived  from  the  principle  of  virtual  work 
given  in  equations  (1)  and  (2).  However,  in  contrast  to  the  total  Lagrangian  formulation,  the 
updated  Lagrangian  method  involves  the  principle  of  virtual  work  being  established  for  the 
equilibrium  configuration  at  the  previous  load  increment  at  time  t. 

(  “7  T’n+l  ^^n+l  dfl  ^n+1  (1^) 

JOt 
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where  Tn+i  and  £„+!  denote  the  second  Piola-Kirchkoff  stress  tensor  and  the  Green-Lagrange 
strain  tensor  measured  in  the  updated  configuration  at  time  t,  respectively.  The  equilibrium 
configuration  at  time  t  is  then  called  the  reference  configuration  for  current  load  step.  The 
updated  Lagrangian  formulation  can  therefore  be  visualized  as  a  piecewise  total  Lagrangian 
formulation. 

In  the  updated  Lagrangian  formulation,  an  incremental  displacement  is  defined  with  re¬ 
spect  to  the  configuration  at  time  t.  The  stress  components  can  therefore  be  expressed  as  the 
sum  of  stress  components  at  time  t  and  a  small  increment  in  the  same  way  as  in  (12).  The 
strain  components  are  similarly  decomposed  into  their  finear  and  nonlinear  parts, 

£n+l  =  ^n+l  +  Vn+l  (20) 

where  the  linear  part  is  defined  as 

6n+l  =  2  ^i,i)n+l 

and  the  non-linear  part  is  defined  as 

‘Hn+l  — 

Remark:  It  is  important  to  note  that  the  initial  displacement  effect  disappears  in  the 

updated  Lagrangian  formulation.  Actually,  it  has  been  taken  into  account  implicitly  through 
the  continuous  updating  of  the  reference  configuration. 

Substituting  the  updated  stresses  and  strains  in  (19),  and  introducing  the  linearized  in¬ 
cremental  constitutive  relation 

cr  =  jT  =  c:e  (21) 

and  neglecting  the  quadratic  terms,  we  obtain  the  linearized  equation  of  motion 


e:c:Sedn,+  tr  •  Sri  dQ  — 


aSedO. 


(22) 


Defining  the  out-of-balance  virtual  work  as 


Residual  =  =  /' 


pext 

^n+l 


/ 


(i) 

n+1 


,  dn 


(23) 


We  establish  the  Newton- Raphson  iteration  procedure  for  updated  Lagrangian  formulation 
as  follows, 


c:  Se^+i  d^l  + 


0*714-1  67] 


(i+1) 

ti4-1 


/ext 

n4-l 


Se%,dn 


(24) 


After  each  iteration,  displacements  are  accumulated  in  the  same  way  as  in  the  total 
Lagrangian  formulation  (7)  and  the  incremental  second  Piola-Kirchhoff  stress  and  Green- 
Lagrangian  strain  tensors  are  evaluated  which  are  then  pushed  forward  to  the  current  con¬ 
figuration.  In  order  to  evaluate  the  out-of-balance  virtual  work  in  (24),  a  transformation  has 
to  be  performed  to  obtain  the  current  Cauchy  stress  tensor  as  follows. 


1 


^714-1 


Jn 


Tn  + 


1 


{Fn+l^Sn+l 


(25) 


The  variation  of  strain  tensor  are  now  evaluated  with  respect  to  the  latest  config¬ 

uration  as 

=  <5^  (Vu(*)  -h  (26) 

and  furthermore,  the  integration  of  the  right  hand  side  of  (24)  is  carried  out  over  the  latest 
configuration  as  well. 

Substituting  now  the  element  geometry  and  displacement  interpolations  into  equations 
(16)  and  (24)  as  we  would  do  in  linear  analysis,  we  obtain  the  incremental-iterative  procedures 
for  general  finite  element  nonlinear  analysis.  For  total  Lagrangian  formulation,  we  have 


Ipext 

^n4-l 


(27) 


and  for  updated  formulation,  we  have, 

(fci. + fc^)A«<‘+')  =  n’ii  -  /sf  (28) 
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In  both  (27)  and  (28),  the  iteration  matrices  called  the  tangent  stiffness  matrices  are  on 
the  left-hand  side  and  include  the  elastic  stiffness  matrix  and  geometric  stiffness  matrix.  The 
right  hand  side  vectors  are  out-of-balance  force  vectors  which  represent  the  difference  between 
the  externally  applied  load  vector  and  the  internal  force  vector  associated  with  the  stresses 
in  the  structure. 

Remark:  It  can  be  shown  that  the  tangent  stiffness  matrices  and  out-of-balance  force  vectors 
evaluated  in  equations  (27)  and  (28)  are  identical.  (See  Bathe  [1994]). 

To  appreciate  the  difficulties  of  implementing  the  Lagrangian  formulations  into  an  existing 
standard  linear  finite  element  code,  it  is  helpful  to  compare  (27)  and  (28)  with  standard  finite 
element  equations  for  linear  analysis,  say 


Kd  =  F 


It  is  readily  recognized  that  in  the  Total  Lagrangian  formulation,  although  all  the  spatial 
derivatives  and  integrations  are  with  respect  to  the  initial  undeformed  configruation  and  the 
geometric  stiffness  matrix  has  the  same  form  as  that  in  linear  buckling  analysis,  we  still  have 
three  major  difficulties: 

(1)  Because  of  the  initial  displacement  effect,  the  formulation  of  the  elastic  stiffness  matrix 
Kli  becomes  much  more  complicated  than  that  in  the  linear  analysis. 

(2)  The  addition  of  the  quadratic  displacement  derivative  terms  in  the  expression  of  the  Green- 
Lagrange  strain  tensor  to  the  linear  strain-displacement  relation  may  require  significant 
changes  in  organization  of  the  linear  finite  element  code. 

(3)  The  evaluation  of  internal  force  vectors  is  computationally  complicated  and  involves  con¬ 
siderable  coding  efforts. 

In  the  updated  Lagrangian  formulation,  because  the  continuous  updating  of  the  reference 
configuration  eliminates  the  initial  displacement  effect,  both  the  elastic  stiffness  matrix, 
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and  the  geometric  stiffness  matrix  have  the  same  form  as  that  in  linear  analysis.  Fur¬ 
thermore,  the  formulation  of  internal  force  vector  is  somewhat  simplified  as  well.  However, 
we  still  have  a  number  of  difficulties  as  enumerated  below: 

(1)  We  have  to  deal  with  quadratic  terms  in  the  Green-Lagrange  strain  tensor. 

(2)  We  need  to  carry  out  a  complicated  transformation  to  establish  the  current  Cauchy  stress 
tensor. 

(3)  We  need  to  calculate  the  new  spatial  derivatives  with  respect  to  the  latest  obtained  con¬ 
figurations  at  the  end  of  each  iteration. 

If  we  want  to  solve  more  general  nonlinear  problems,  such  as  those  involving  large  strains  or 
nonlinear  materials,  the  exact  Lagrangian  approaches  must  be  employed  to  generate  reliable 
finite  element  solutions.  Fortunately,  in  the  case  of  small  strains,  a  series  of  approximations 
can  still  be  made  to  degenerate  the  updated  Lagrangian  formulation  to  element  co-rotational 
procedure,  as  explained  below. 

2.3  Element  Corotational  Procedure 

The  total  and  updated  Lagrangian  formulations  discussed  in  previous  section  are  so  general 
that  they  are  able  to  handle  any  kinds  of  nonlinearity,  i.e.,  large  displacements,  large  rotations, 
finite  strains  as  well  as  various  forms  of  nonlinear  constitutive  relations.  However  in  the  case  of 
small  strains,  a  series  of  approximations  can  be  made  to  significantly  simphfy  the  Lagrangian 
formulations. 

Theoretically,  an  arbitrary  motion  of  a  general  continuous  medium  can  always  be  decom¬ 
posed  into  a  rigid  body  motion  followed  by  a  pure  relative  deformation.  In  finite  element 
analysis,  this  decomposition  can  be  accomplished  by  defining  a  local  convected  coordinate 
system  for  each  element  which  translates  and  rotates  with  the  element,  but  not  deform  with 
the  element.  If  the  finite  element  is  sufficiently  small,  the  pure  deformation  part  of  displace¬ 
ment  obtained  by  subtracting  the  rigid-body  motion  components  from  the  total  displacement 
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is  normally  small.  Consequently,  the  linear  theory  can  be  apphed.  This  is  the  key  idea  of  the 
“Co-rotational  Procedure” . 

In  the  following,  we  endeavour  to  show  the  close  relationship  between  co-rotational  proce¬ 
dure  and  updated  Lagrangian  formulation.  Prom  this  relationship  we  establish  some  simple 
but  mathematically  consistent  procedures  for  updating  element  stresses  and  calculating  in¬ 
ternal  force  vector. 

Because  of  the  small  strain  assumption,  the  pure  deformation  part  of  the  displacements  is 
always  a  small  quantity  as  compared  to  the  element  dimensions  in  the  local  convected  coordi¬ 
nate  system.  Consequently,  it  is  reasonable  for  us  to  make  the  following  two  approximations: 

1)  The  changes  in  element  shapes  are  small  in  each  individual  local  increment;  and 

2)  The  displacement  derivatives  of  the  pure  deformation  part  of  the  displacements  with 
respect  to  the  current  configuration  measured  in  the  local  convected  coordinate  system 
axe  of  the  same  order  as  the  small  strains.  Furthermore,  the  quadratic  terms  in  the 
Green-Lagrange  strain  components  are  neglected. 

The  second  approximation  was  first  introduced  by  Belytschko  and  Hsieh  [4].  Prom  the 
above  two  approximations,  it  follows  that  the  linear  theory  can  be  employed  in  the  element 
local  convected  coordinate  system  to  provide  geometrically  nonlinear  solutions  to  large  dis¬ 
placements,  large  rotation,  but  small  strain  problems. 

Remark:  Although  all  the  motions  of  a  continuous  body  can  be  theoretically  decomposed 
into  rigid-body  and  pure  deformation  components,  in  practical  computations,  the  numerical 
performance  of  the  co-rotational  procedure  greatly  depends  upon  the  method  of  construction 
of  the  local  convected  coordinate  system.  The  more  completely  the  rigid  body  motion  is  iso¬ 
lated  from  the  total  displacement  field,  the  more  accurate  the  finite  element  solution  provided 
by  the  co-rotational  procedure  is. 

In  order  to  present  the  corotational  procedure,  we  refer  the  three  sets  of  coordinate  sys- 
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terns.  They  are: 

(i)  The  Global  Coordinate  System  which  is  the  Cartesian  coordinate  system  fixed  in  space. 

(ii)  The  Local  Convected  Coordinate  System,  I  =  [ii,i2>*3]  which  is  associated  with  each 
element  and  undergoes  rigid  body  rotations  and  translations  with  the  element,  and 

(iii)  The  Surface  Coordinate  System,  Vi,  which  is  the  nodal  triad  associated  with  node  i  in  a 
structural  element.  During  the  equilibrium  iteration,  Vi  experiences  both  the  rigid  body 
rotations  and  pure  deformation.  Obviously  Vi  can  be  measured  in  either  global  or  local 
convected  coordinate  system,  as 

Vi  =  =  I  =  [Vi,l,i>i,2,Vi,3] 

where  i  denotes  the  node  number  and  tilde  indicates  quantities  measured  in  the  local 
convected  coordinate  systems. 

In  order  for  the  linear  finite  element  theory  to  generate  sufficiently  accurate  results  for 
nonlinear  problems,  the  deformation  part  of  the  displ2icements  must  be  kept  as  small  as 
possible  during  the  solution  process.  As  a  result,  it  is  appropriate  to  update  the  reference 
configuration  after  the  equilibrium  state  has  been  reached  for  each  load  increment.  In  this 
way,  the  displacement  we  are  dealing  with  will  be  just  the  displacement  increment  in  the 
current  load  step  defined  in  equation  (17).  Because  this  concept  is  identical  to  that  employed 
in  the  updated  Lagrangian  formulation,  equation  (28)  wiU  be  an  appropriate  starting  point 
for  our  derivation  of  the  co-rotational  procedure. 

It  has  been  pointed  out  in  the  previous  section  that  in  updated  Lagrangian  formulation,  if 
the  modified  Newton-Raphson  iterative  method  is  adopted,  the  elastic  and  geometric  stiffness 
matrices  have  the  same  form  as  in  linear  analysis.  (See  Bathe  [1994]).  The  main  difficulty 
arises  in  the  calculation  of  the  internal  force  vector.  The  internal  force  vector  in  (28)  can  be 
expressed,  for  an  individual  element,  as: 
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where  has  the  same  form  as  in  ordinary  hnear  finite  element  analysis,  but  contains  the 

spatial  derivatives  with  respect  to  the  latest  configuration  obtained,  i.e.,  which  is  not 

necessarily  in  equihbrium  state.  is  the  Cauchy  stress  tensor  and  is  also  measured  in 

the  latest  configuration.  To  obtain  exact  Cauchy  stress  tensor,  we  have  to  first  calculate  the 
incremental  Green-Lagrange  strain  from  which  the  incremental  second  Piola-Kirchhoff  stress 
tensor  can  be  obtained,  and  secondly,  perform  the  complex  stress  transformation  as  described 
in  proceeding  section.  All  variables  in  (29)  are  with  respect  to  the  global  coordinate  system. 

To  apply  the  corotational  procedure,  it  is  convenient  to  transform  equation  (29)  into  local 
convected  coordinate  system,  as 

(30) 

•'^n+l 

where  the  tilde  indicates  quantities  measrued  in  the  local  convected  coordinate  system. 
represents  the  strain-displacement  matrix  in  the  local  convected  coordinate  system  in  config¬ 
uration 

The  Cauchy  stress  tensor  in  the  latest  obtained  configuration  can  be  calculated  through 
the  stress  transformation  in  local  convected  coordinate  system  as  follows 

Jo-  =  f  =  FSF’^  (31) 

a  =  ^FSF'^  (32) 

tJ 

Consequently,  stress  update  in  the  current  configuration  can  be  expressed  as 

ff„+i  =  ffn  +  if'ASF’’  (33) 

Based  on  the  first  assumption  where  the  changes  on  element  shapes  are  assumed  to  be 
negligible,  the  deformation  gradient  matrix  in  local  convected  coordinate  system  must  be 


close  to  an  identity  matrix  and  the  change  in  mass  density  is  also  negligible,  we  have  Po  =  Pt 
for  J  =  1  and  the  deformation  gradient 

F.  -  1  ^ 

~  ”  I*  0  ^  ^ 

Substituting  (34)  into  (33),  we  get  further  approximation 

d^n+l  ~ 

=  +  C: 

where  because  of  the  small  strain  assumption  in  the  co-rotational  framework,  we  drop  the 
nonlinear  term  in  the  Green  Lagrange  strain  tensor  and  are  left  with  the  linear  strain  tensor. 
Consequently, 


(34) 


(35) 


=  +  C:  i(Vti  +  V*’’)<‘i, 

=  ?(<>  +  C:  (36) 


where  we  have  used  the  transformation  relation  between  the  nodal  degrees  of  freedom  in  the 
global  and  the  local  convected  coordinate  systems.  Once  the  internal  force  vector  is  obtained 
in  local  convected  coordinate  system,  a  simple  coordinate  transformation  can  be  performed 
to  obtain  the  internal  force  vector  in  global  coordinate  system  as: 


_  jT 

Jn+1  —  ■*  Jn+1 


(37) 


where  the  transformation  matrix  is  obtained  as 


mf 

n\ 

l\mi 

mini 

nih 

q 

ml 

nl 

l2m2 

m2n2 

n2^2 

^3 

ml 

nl 

msna 

nsh 

2I1I2 

2m\m2 

2nin2 

l\m2  +  ^2^1 

min2  +  m2ni 

n\l2  +  n2ii 

2/2^3 

2m2Tn3 

2n2Ti3 

^2^3  d"  ^3UT’2 

m2n3  +  m^n2 

n2h  +  nzl2 

2hh 

2m^mi 

2n3ni 

h'^i  +  hmz 

m^rii  +  miUs 

n^h  +  ni/3 

(38) 
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and 


11  =cos(i,ei)  mi  =  cos(jf,ei)  ni  =cos(fc,ei) 

12  =  cos{i,  €2)  m2  =  cos(j,  62)  n2  =  cos(fc,  62) 

Z3  =  cos(i,  €3)  m3  =  cos(j,  63)  n3  =  cos(fc,  63) 

The  relation  between  the  Cauchy  stress  tensor  in  the  global  and  the  local  convected 
coordinate  system  is  expressed  through  the  standard  transformation  matrix  as 

^i+i  =  (39) 


Remark:  The  above  treatment  to  the  geometric  nonlinearites  is  element-independent  and 

applicable  to  all  finite  elements.  However,  the  definition  of  the  local  convected  coordinate 
system  is  still  element-dependent  and  special  care  has  to  be  excercised  to  obtain  a  better 
numerical  performance.  Another  element-  dependent  step  is  the  modification  of  nodal  tri¬ 
ads  after  each  iteration.  The  formulations  must  be  consistent  with  the  kinematic  relations 
employed  in  the  element. 
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3.  STABILIZATION  TECHNIQUES  FOR  UNDERINTEGRATED 


NONLINEAR  SHELL  ELEMENTS 
3.1  Introduction 

In  large-scale  finite  element  analyses,  a  significant  number  of  elements  and  large  computer 
memory  are  required  to  obtain  the  detailed  information  for  engineering  design  or  process 
control.  When  these  analyses  are  explicit,  computational  costs  axe  mostly  determined  by  the 
efficiency  of  the  elements,  especially  for  non-linear  problems.  In  the  past,  full  quadrature 
and  selective  reduced  integration,  were  extensively  used  in  the  finite  element  analyses.  Full 
quadrature  ensures  the  stability  and  convergence  of  solutions.  However,  it  is  very  expen¬ 
sive  owing  to  the  requirement  of  many  computational  operations  for  evaluating  the  element 
stiffness  matrices  and  the  internal  force  vectors  and  it  also  lead  to  volumetric  locking  and 
transverse  shear  locking  in  thin  structure  bending  and  incompressible  problems.  A  remedy 
for  this  is  to  use  the  selective  reduced  integration  in  which  the  full  quadrature  and  reduced 
quadrature  are  applied  to  different  terms  to  form  the  element  stiffness  matrix  as  proposed  by 
Malkus  and  Hughes  [14]  and  Nagtegaal  et  al.  [15],  among  others.  However,  it  is  as  costly  as 
full  quadrature. 

Perhaps,  the  most  efficient  elements  are  the  one-point-quadrature  elements  with  hourglass 
control  developed  by  Flanagan  and  Belytschko  [16],  Belytschko  [17]  and  Belytschko  et  al.  [18]. 
The  mesh  instability  associated  with  the  underintegrated  elements  is  controlled  by  adding  a 
stabilization  to  the  one-point-quadrature  element.  The  stabilization  terms  are  developed  to 
ensure  the  consistency  of  the  finite  element  equations  and  its  magnitude  is  controlled  by 
a  user  input  stabilization  parameter.  Liu  and  Belytschko  [19]  also  developed  a  one-point- 
quadrature  element  for  heat  conduction  problems.  In  this  work,  the  stabilization  parameter 
is  determined  by  solving  an  eigenvalue  problem.  In  Beljdschko  et  al.  [20]  the  Hu-Washizu 
variational  principle  is  used  to  examine  the  magnitude  of  the  stabihzation  parameters.  More 
recently,  an  assumed  strain  stabihzation  of  the  four-node  quadrilateral  element  and  the  eight- 
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node  hexahedral  element  with  one-point  quadrature,  where  the  stabihzation  parameters  are 
not  required,  was  proposed  by  Belytschko  and  Bindeman  [20,21]. 

An  alternative  approach  for  hourglass  control  is  proposed  by  Liu  et  al.  [22]  in  which  the 
resulting  stabilization  matrix  requires  no  stabilization  parameter.  It  is  shown  that  the  stabi¬ 
lization  vector  7  can  be  obtained  simply  by  taking  the  partial  derivatives  of  the  generalized 
strain  vector  with  respect  to  the  natmal  coordinates.  The  strain  vector  is  therefore  approx¬ 
imated  by  the  combination  of  a  constant  part  and  other  parts  involving  strain  derivatives. 
However,  shear-related  locking  phenomena  are  not  taken  into  consideration  and  no  three  di¬ 
mensional  result  is  reported  in  their  study.  Belytschko  [18]  established  an  hourglass  control 
procedure  by  expanding  the  stress  in  a  Taylor  series  about  the  element  center.  It  is  found 
that  the  hourglass  control  is  suppressed  by  the  retention  of  the  first  and  second  derivatives 
terms  in  the  expansion.  The  detailed  derivation  of  element  matrices  and  numerical  examples 
are,  however,  only  for  two  dimensional  problems. 

In  this  work  a  three-dimensional  underintegrated  element  based  on  the  procedures  similar 
to  that  proposed  by  Liu  et  al  [22]  are  developed.  Emphasis  is  placed  on  avoiding  of  locking 
and  the  removal  of  spurious  singular  modes.  This  element  is  applicable  to  plate  and  shell 
bending  problems,  and  more  importantly  it  is  suitable  for  extension  to  elasto-plasticity.  Sec¬ 
tion  3.2  presents  the  synopsis  of  the  underintegrated  four- node  quadrilateral  elements  with 
hourglass  control  as  developed  by  Belytschko  et  al .  The  three-dimensional  element  suitable 
to  bulk  deformation  and  plate/shell  structural  analyses  is  developed  in  Section  3.3. 
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3.2  Four-node  quadrilateral  element 


Fig.  1.  Physical  and  reference  configurations. 


Let  us  consider  a  four- node  quadrilateral  element  as  shown  in  Figure  1.  The  configuration 
of  the  element  is  a  bi-unit  square  in  the  natural  coordinate  system  The  spatial  coordi¬ 

nates  X  and  the  velocity  component  v  in  the  element  are  approximated  by  linear  combinations 
of  nodal  values  and  and  the  element  shape  functions. 

The  gradient  submatrix  is  given  by 


Ba 


^a,x 

N 

^ya,y 


which  can  be  evaluated  at  the  element  center  as 


Ba(0) 


>a.x  (0)' 

bia 

(0) 

b2a 

where 


^>1  =  {^>io}  =  ^[2/24  »  ysi ,  2/42  )  yis]^ 


62  ==  {^20}  =  ^[^42  ,  2^13  )  2:24  ,  2:31 


(40) 


(41) 


(42) 

(43) 


XIJ  =  Xj  -Xj 
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(44) 


Fig.  2.  Bending  Modes. 


yij  =  yi  -yj  (45) 

=  2(^31  y42  +  ^24y3l)  (46) 

and  /  and  J  designate  the  node  numbers. 

For  the  purpose  of  identifying  the  deformation  modes  of  element,  let  us  define  column 
vectors 


s  =  [1,1,1,  iF 

(47) 

=  Xi  =  [xi  ,  X2  ,  X3  ,  3:4]^ 

(48) 

=  X2  =  [yi ,  y2  ,  2/3 , 2/4]^ 

(49) 

h  =  [1,-1, 1,-1]^’ 

(50) 

€  =  [-1, 1, 1,  -iF 

(51) 

ri  =  [-1,-1,  l,lF 

(52) 

where  x  and  y  are  the  element  nodal  coordinates  of  the  physical  space;  ^  and  rj  are  element 
nodal  coordinates  of  the  biunit  square:  and  h  is  the  hourglass  vector.  In  Figure  2,  we  show 
two  bending  modes:  {h,  0}  and  {0,  h}  which  are  the  deformation  mode  associated  with  no 
energy  in  the  one-point  quadrature  element  but  resulting  in  a  non-constant  strain  in  the 
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element. 


The  strain  rate  i  is  approximated  by  expanding  it  in  a  Taylor  series  about  the  element 


center  up  to  linear  terms: 


=  HO)  +  e,j(0)f  +  £,,(0)i) 


(53) 


or 


where 


rien 

e  =  Ba{^,'n)Va 

a=l 


=  Ba{0)  +  Ba,^{0)^  +  Ba,r,{0)v 


(54) 


(55) 


The  first  term  on  the  right  hand  side  of  equation  (53)  is  the  constant  strain  rate  evaluated  at 
the  quadrature  point,  0,  and  the  other  terms  are  linear  strain  rate  terms. 

After  some  algebra,  the  explicit  expressions  of  the  first  derivatives  of  Ba{0)  with  respect 
to  natural  coordinates,  Ba,^{0)  and  Ba,r]{0),  can  be  shown  to  be 


where 


where 


(0) 


bl,^a 

52, {a 


Ba,r,  (0)  = 


51, jja 

52, Tja 


bi,^  =  =0127 

b2,^  =  {52, =oii7 
5i,»j  =  {5i,,,J  =  0227 
52,7/  =  {52,T)a}  =  “2l7 


Oil 

Ol2 


IT 

-ey 

4A 


(56) 

(57) 

(58) 

(59) 

(60) 
(61) 

(62) 

(63) 
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(64) 


0-21 


—rj^x 

4A 


The  7  in  equations  (58)-(65) 
of  B(0).  It  is  given  by 


a22  = 


■n'^y 

AA 


(65) 


is  the  stabilization  vector  which  spans  the  improper  null-space 


7  =  /i  -  ( h'^Xi)  bi 


(66) 


where  i  is  the  summation  index  from  1  to  2.  7  is  orthogonal  to  the  linear  velocity  field  and 
provides  the  proper  stabilization  for  the  element.  Belytschko  and  co-workers  [17,20]  derived 
these  vectors  from  the  requirements  of  consistency  of  the  finite  element  equations  in  the  sense 
that  the  gradients  of  the  linear  fields  were  evaluated  correctly. 

To  alleviate  volumetric  locking,  we  employ  the  ideas  underlying  selective/reduced  inte¬ 
gration.  Ba{^,r])  is  decomposed  into  two  parts  as 


S.(e.*?)  =  -B"'(0)  +  Srte>))  (67) 

where  are  the  gradient  submatrices  due  to  the  dilatational  part  of  Ba,  and  are 

the  gradient  submatrices  due  to  the  deviatoric  part  of  Ba  ■  Here,  the  dilatational  part  of  the 
gradient  matrix  has  been  underintegrated  and  evaluated  only  at  one  quadrature  point,  0,  to 
alleviate  the  volumetric  locking.  Expanding  about  the  element  center  via  equation  (55), 
equation  (67)  can  be  written  as 

Bad.ri)  =  B.(0)  +  B^S'(O)^  +  (68) 

where  Ba(0)  are  the  one-point-quadrature  gradient  submatrices  contributed  from  both  the 
dilatational  and  deviatoric  parts.  The  other  terms  on  the  right-hand  side  of  above  equation 
are  the  gradient  submatrices  corresponding  to  non-constant  deviatoric  strain  rates.  It  is  noted 
that  the  resulting  element  using  the  gradient  matrices  as  in  equation  (67)  or  (68)  exhibits 
no  hourglass  mode  if  the  element  internal  virtual  work  is  evaluated  by  using  a  multiple- 
quadratvue-point  integration. 
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The  element  developed  so  far  is  not  suitable  for  the  thin  beam  or  plate  analysis  owing 
to  the  so  called  shear  locking  for  the  bending  modes.  It  is  well  known  that  the  reduced 
integration  on  shear  strain  rates  in  the  global  coordinate  system  is  not  physically  correct 
when  the  referential  coordinate  system  is  not  aligned  with  the  global  one.  Hence,  a  co- 
rotational  coordinate  system  which  rotates  with  the  element  is  used  for  the  derivation  that 
follows.  To  alleviate  the  shear  locking,  the  gradient  submatrices  for  the  general  2-D  element 
are  interpolated  in  the  following  form: 


B..  («,>?)  =  (0)  +  (0)  I  +  sj-  (0)  n  (69) 

Byy  it  v)  =  B,y  (0)  +  (0)  (  +  (0)  ,  (70) 

({.»))  =  Bx„(0)  (71) 

B„  (?. ,)  =  Bf-j  (0)  {  +  B.^"  (0)  V  (72) 

where 


Bxx  (0) 

bT 

0 

Byy  (0) 

0 

b^ 

Bxy  (0) 

bT 

B,,  (0)^ 

0 

0 

Bj-J  (0) 

§012  7^ 

-5^117^ 

(») 

-iai2  7^ 

fan  7^ 

BS  (0) 

-iai2  7^ 

1  T 

B^y  (0) 

i«22  7'^ 

-|a2i7’’ 

c,  (“) 

= 

1  T 

-5^22  7 

ia2i7^ 

bSj,  (0) 

1  T 

-3^22  7 

-|a2i7^ 

Bxx  5  Byy  ,  Bxy  and  are  the  gradient  submatrices  corresponding  to  the  shear  strain  rates 
ixx,£yy,£xy  and  izz,  respectively.  Here,  only  the  constant  term  is  used  for  the  shear  strain 
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rate  component  such  that  the  mode  causing  locking  is  removed.  The  normal  strain  rates  keep 
all  non-constant  terms  given  in  equation  (68). 

Although  one-point  quadrature  with  hourglass  control  can  be  applied  to  the  element 
developed  herein,  it  might  be  not  accurate  enough  for  more  advanced  nonhnear  constitutive 
models.  Hence  we  follow  Liu  et  al.  [22]  to  use  two-point  quadrature  where  the  element 
internal  force  vector  is  evaluated  at  the  two  integration  points  located  as  follows: 

Pointl  . 

Pomt2  : 

This  two-point  quadrature  element  exhibits  no  hourglass  mode  and  is  rank  sufficient.  By 
assuming  that  the  Jacobian  is  a  constant,  one  half  of  the  area,  the  element  internal  force 
vector  can  be  integrated  as  follows: 

/“  =  E  I®’’ K*) 

1=1 

where  denotes  the  natural  coordinates  of  the  integration  point  i  and  A  is  the  element  area. 
The  element  internal  force  vector  can  be  further  rearranged  as 


(77) 


where  and  internal  force  vector  resulting  from  one-point  quadrature  and 

stabilization,  respectively.  They  are  given  by 


and 


j:int  _  A  fiibi  +  fi2b2 

^  2  '5^12^1  +  T22  ^2 

(ai2  +  022)  X  (2tii  —  7)22  —  ^33)  7 
(oii  +  021)  X  (-fii  +  2f22  -  fas)  7 


(78) 

(79) 
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where  aij  are  given  in  equations  (62)-(65)  and 


=  Tij  (6)  +  Tij  (42) 

(80) 

=  Tij  (6)  -  rij  (^2) 

(81) 

3.3  Eight-node  Hexahedral  element 


Fig.  3.  Physical  and  reference  configurations  for  8-node  element. 


Let  us  consider  an  eight-node  hexahedral  element  as  shown  in  Figure  3.  The  spatial 
coordinates  and  velocity  in  the  element  are  approximated  by  linear  combination  of  nodal 
values  and  shape  functions  as  follows: 


ricn 


^  ^  (^en  —  8) 

a=l 

(82) 

Vi  =  ^  Na{^,T],0Via 

(83) 

a=l 


C)  —  g  ^“0  (1  +  VaV)  (1  +  CaC) 
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(84) 


where  the  subscripts  i  and  a  denote  coordinate  component  ranging  from  one  to  three  and  the 
element  node  number,  ranging  from  one  to  eight,  respectively. 

For  the  purpose  of  identifying  the  deformation  modes  of  the  element,  let  us  define  the 
gradient  submatrices  Ba{0)  and  other  column  vectors  as 


(0) 

bla 

Ba  (0)  = 

N..,  (0) 

b2a 

JVa,.  (0) 

&3a 

(85) 


=  [1,1, 1,1, 1,1, 1,1]  (86) 

=  [xi,X2,X3,X4,X5,Xe,X7,X8]  (87) 

=  bi,y2,z/3,2/4,y5,y6,y7,2/8]  (88) 

z'^  =  [ZI,Z2,Z3,Z4,Z5,Z6,Z7,Z8\  (89) 

hj  =  [1,-1, 1,-1, 1,-1, 1,-1]  (90) 

h'^  =  [1,-1, -1,1, -1,1, 1,-1]  .  (91) 

hj  =  [1,1, -1,-1, -1,-1, 1,1]  (92) 

hj  =  [-1,1, -1,1, 1,-1, 1,-1]  (93) 

=  [-1,1, 1,-1, -1,1, 1,-1]  (94) 

=  [-1,-1, 1,1, -1,-1, 1,1]  (95) 

=  [-1,-1, -1,-1, 1,1, 1,1]  (96) 


where  x,  y  and  z  are  the  nodal  coordinates  and  hi  is  the  ^77-hourglass  vector,  /12  the 
hourglass  vector,  /13  the  ryC-hourglass  vector  and  /14  the  ^77C-hourglass  vector. 


The  Jacobian  matrix  evaluated  in  the  center  of  an  element  can  be  shown  to  be 

1 


J(0)  =  [Ji^]  = 


8 


$.'^x  i'^z 

rf'x  rj^y  rf^z 

C'^x  Cy  C'^z 


hJ  =  1>2,3 


(97) 
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The  determinant  of  the  Jacobian  matrix  is  denoted  by  Jq  and  the  inverse  of  J(0)  is  given  by 

D: 

D  =  lAjl  =  J-'  (0)  (98) 

The  gradient  vectors  61,62  aiid  63  (which  are  evaluated  at  0)  in  equation  (85),  can  be 
shown  to  be 

61  =  {6i„}  =  +  T)i277  +  r>i3C]  (99) 

62  =  {62a}  =  ^p2i^  +  D22ri  +  T>23C]  (100) 

63  =  {^3a}  =  g  [-031^  +  D^2'n  +  ■^^33C]  (101) 

Expanding  strain  rate  e  in  Taylor  series  about  the  element  center 

£(^)^jC)  =  ^(0)  +  £,{(0)^  +  +  ^,c(o)C 

+  (0)  (0)  (0)  (102) 

or 

rien 

^  =  '^Ba{i,r},C)Va  (103) 

a=l 

where 

=  -B^O)  +  B.,((0)e  +  B„(0),  +  B.,ao)< 

+  B.,j,  (0)  if,  +  (0)  i)C  +  B,,;j  (0) «  (104) 

The  first  term  on  the  right  hand  side  of  the  equation  (102)  or  (104)  is  the  constant 
strain  rates  evaluated  at  the  quadrature  point,  0,  and  the  remaining  terms  are  the  linear  and 
bilinear  strain  rate  terms.  After  some  tedious  algebra,  it  can  be  shown  that  the  first  and 
second  derivatives  of  .Ba(O)  with  respect  to  natural  coordinates  are  given  by 

=  {^a,x^}  =  g  [.Di27i  +  £>1372]  (105) 

62,$  =  {Na^}  =  g  [-D2271  +  -^>2372]  (106) 
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where 


=  {Na,zi}  =  g  [-03271  +  -03372]  (107) 

bl,r)  =  {Na,xr)}  =  g  [Oii7i  +  O1373]  (108) 

^>2,7?  =  {^a,yT)}  =  g  [O2171  +  O2373]  (109) 

^3,7J  =  {-^a,zrj}  =  g  [O3171  +  O3373]  (110) 

bl,C  =  {-^a.ic}  =  g  [Oll72  +  O1273]  (111) 

^>2,C  =  {^a,yc}  =  g  [O2172  +  O2273]  (112) 

bsx  =  {^a,zc]  =  g  [O3172  +  O3273]  (113) 

Kin  =  {^aMn}  =  ^[Oi3  74  -  ip1xi)bi^^  -  (rfxi)6i,^]  (114) 

b2,in  =  {Na,yir,}  =  ^[023  74  “  (P2  ^^i)  “  (^^a^x)  (US) 

b3,ir,  =  {Na,zin}  =  ^[033  74  “  {pI  ^i)  -  (r^JCi)  6i,„]  (116) 

bi^nC  —  {-^a,TJ7c}  ~  g  [Oil  74  ~  Ki  ^j)bi^r}  ~  (Pi  aJx)  hx,^]  (117) 

b2,r]C,  ~  {-^aij/Tic}  ~  g  [O21 74  ~  (^2  ^i)  bi,r]  ~  {P2^i)bix]  (113) 

b3,nC  =  {Na,zr,c}  =  ^[03174  -  {q3^i)bi,r,  “  (Ps  ^i.cl  (HO) 

hl.^C  ~  {^a,xic}  —  g[Oi2  74  ~  {Ql^i)biX  ~  (^1  aix)  ^i.cl  (120) 

b2XC  =  {^a,yic}  =  ^[022  74  “  (ql’  ®x)  bx,^  "  (^’I’a:,)  bx,c]  (121) 

b3,ii  =  {Na.zic)  =  ^[032  74  "  {ql  ^i)  ”  (T’l’aJx)  6x,c]  (122) 

Pi  =  Dahl  +  Dishs  (123) 

qi  =  Oxi/i2  +  Ox2h.3  (124) 

Tj  =  Di2hi  +  Oi3/l2  (125) 
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The  7a  in  equations  (105)-(122)  are  the  stabilization  vectors  which  span  the  improper  null- 
space  of  B(0).  They  are  given  by 

7a  ^  ha  -  {h'^Xi)bi  (126) 

where  i  is  the  summation  index  from  1  to  3.  7a  is  orthogonal  to  the  hnear  displacement  field 
and  provide  the  proper  stabilization  for  the  element. 

Similar  to  the  four-node  quadrilateral  element,  Ba  (^,  C)  i®  decomposed  into  two  parts, 
dilatational  part  and  the  deviatoric  part.  The  dilatational  part  of  gradient  matrices  are 
under  integrated  and  evaluated  at  only  one  quadrature  point,  0,  to  avoid  volumetric  locking: 

((,  n,  C)  =  Sj"  (0)  +  "  tt,  n,  0  (127) 

Expanding  Bf^''  about  the  element  center,  equation  (126)  can  be  written  as 

Ba  (^ ,  n,  C)  =  Ba  (0)  +  B„7  (0)  C  +  (0)  r,  +  (0)  C 

+  BJs;  (0)  in  +  B^;^  (o)  n(  +  (o) «  (i28) 

where  Ba  (0)  are  the  one-point-quadrature  gradient  submatrices  contributed  from  both  the 
dilatational  and  deviatoric  parts.  The  remaining  terms  on  the  right-hand-side  of  the  above 
equation  are  the  gradient  submatrices  corresponding  to  non-constant  deviatoric  strain  rates. 
It  is  noted  that  the  element  using  the  gradient  matrices  as  in  equation  (127)  or  (128)  is 
properly  underintegrated  and  exhibits  no  hourglass  mode  if  the  element  internal  energy  is 
evaluated  by  using  the  multiple-point  quadrature. 

The  element  developed  so  far  is  not  suitable  to  plate  /shell  analysis  owing  to  the  shear 
and  membrane  locking  in  thin  structures.  To  remove  shear  locking,  the  gradient  submatrices, 
corresponding  to  the  assumed  shear  strain  rates  is  written  in  an  orthogonal  co-rotational 
coordinate  system  rotating  with  the  element  as 

Bxx«,>),0  =  Baa(o)  +  bS^(o)?  +  Bj;y,(o)>)  +  Bj;yj(o)c 

+  bS^j,  (0)  in  +  (0)  <  +  Bf-„  (0)  a 
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(129) 


(130) 


where 


^yy  (^)  Vj  C) 


(0)  +  B-";  (0) «  +  B^-  (0) ,  +  Bjjy^  (0)  C 

+  (0)  (0)  IC  +  SJS:«  (0)  « 


B„(f,ri,C)  =  B„(0)  + 


+  Bjyj,  (0)  in  +  Bt:„i  (0)  i,c  +  (0)  ce  (i3i) 


“  Bx„(o)  +  B^;y.(o)c 

(132) 

B,,(f,,),c)  =  B„,  (0)  +  Bj^s  (0)  f 

(133) 

B„(?,i),C)  =  B„(0)  +  B*y,(o), 

(134) 

Bxx  (0) 

bf 

0 

0 

Byy  (0) 

0 

bl 

0 

B,,  (0) 

0 

0 

bl 

Bxy  (0) 

bl 

0 

By,  (0) 

0 

bl 

B„(0) 

1 

COh^ 

0 

b^ 

(135) 


'Bf?f(0)' 

0 

-i^227r 

-1^3372^ 

(0) 

0 

ID221I 

BSi  (») 

1 

0 

-\D221^ 

1^3372"’ 

(0) 

8 

D22'yT 

0 

0 

«ji:£(o) 

0 

£>3372^ 

I?227r 

BSi  (0) 

-D337J 

0 

0 
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BS:,  (0) 

iDnlT  0 

B^.r,  (») 

-{DnlT  0 

-|£»3373"' 

BiT,r,  (0) 

1 

-\DnlJ  0 

(0) 

8 

0  £>ii7i’ 

0 

0  ^337j 

0 

(0) 

£>3371’  0 

£>ii7i’ 

(137) 


'BfzM 

|£>ii7j’ 

-P2273^ 

0 

Bt7.i  (0) 

-|£>n7i’ 

§£>2273^ 

0 

BS  («) 

1 

-5^ii7i^ 

0 

bs:<  (0) 

8 

£>227^ 

Dial 

0 

BtZi  (0) 

0 

0 

£>2273’ 

BJlTc  (0) 

0 

0 

'bs:«,(o)' 

0 

0 

-|i^337r‘ 

(0) 

0 

0 

bs:s,  (0) 

1 

0 

0 

iDss^I 

(0) 

8 

0 

£>337! 

0 

0 

0 

0 

bS,  (0) 

£>3371 

0 

0 

(139) 


’b5"c(o)' 

§£>1171 

0 

0 

Bjs:,c  (0) 

-i£>ii7r 

0 

0 

BS:,c  (0) 

1 

-|£>ii7r 

0 

0 

B^,(  (0) 

8 

0 

£>117J 

0 

B^r,£  (0) 

0 

0 

0 

BS:,c  (0) 

0 

0 

£>117J_ 

(140) 
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(141) 


0 

-\D221I 

0 

s„1:a(0) 

0 

lD22ll 

0 

(0) 

1 

0 

0 

(0) 

8 

022^1 

0 

0 

sS:ce(0) 

0 

0 

D221I 

BiTxi  (0) 

0 

0 

0 

In  above  equations,  B^x,  Byy,  B^z,  B^y,  By^  and  B^x  are  the  gradient  submatrices  cor¬ 
responding  to  strain  rates  ixx,^yy,£zz,^xy,£yz  and  izx  respectively.  Here,  only  one  non¬ 
constant  term  is  used  for  each  shear  strain  rate  component  such  that  the  modes  causing 
shear  locking  are  removed.  The  normal  strain  rates  keep  all  non-constant  terms  given  in 
equation  (127).  In  S‘*®''(0),  only  those  terms  corresponding  to  a  parallelpiped  element  are 
used  for  stabiUzation. 

For  advanced  nonlinear  constitutive  models  in  3-D  calculations  of  large  deformation  prob¬ 
lems,  we  propose  to  use  a  four-point-  quadrature  scheme  instead  of  one-point  quadrature. 
Following  Belytschko  et.  al.  [20]  the  element  internal  force  vector  is  evaluated  at  the  four 
integration  points  located  as  follows; 


Point  1  : 


Point  2  : 


_1_ 

V3’ 


Point  3  ; 


Point  4  ; 


(142) 

(143) 

(144) 

(145) 


By  assuming  that  the  Jacobian  is  a  constant,  one-quarter  of  the  element  voliune,  the 


element  internal  force  vector  can  be  integrated  as  follows: 

/'“  =  E 

i—\ 

where  denotes  the  natural  coordinates  of  the  integration  point  i  and  V  is  the  element 
volume. 

The  element  internal  force  vector  can  be  rearranged  in  the  form 

/“•  =  /i“  +  /iSb  (147) 


where  /J"*  and  are  the  internal  force  vectors  resulting  from  a  four-point  quadrature  and 
the  stabilization  procedure,  respectively.  They  are  given  by 


4  y  +  ri2(Ct)h2  +  Tl3i^i)b3 

~  ^2  X  '^2l(^i)hl  +  T22{ii)b2  +  T23(^t)63 

T3l{^i)bi  +  T32{^i)b2  +  T33(^i)63 


and 


J  stab 


V- 

^  32 


i=l 


+  54(6)'^12(^i)  +  p5(Ci)'^3l(^i) 

92{^i)T22'' i^i)  +  96i^i)Tl2{^i)  +  97{^i)T23{^i) 

93{^i)rii''{^i)  +  &8(6)T23te)  +  59(^i)T3l(6) 


(148) 


(149) 


where  the  superscript,  dev,  denotes  the  deviatoric  part  of  the  stress,  and  the  other  quantities 


are  given  by 


9i  (^)  =  -^11  (’771  +  ^72  +  277C74) 

(150) 

92  (^)  =  D22  (^7i  +  C73  +  2^C74) 

(161) 

93  (€)  =  D33  (^72  +  7773  +  2^7774) 

(152) 

^4(^)  =  D22CI3 

(153) 

93  i^)  =  -D33’?73 

(154) 
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(155) 


97  (^) 

=  -^^33^72 

(156) 

98  (^) 

=  1^22^71 

(157) 

99 

=  DnVlfi 

(158) 

It  is  noted  that  the  element  developed  above  cannot  pass  the  patch  test  because  with  the 
one-point  quadrature  the  element  internal  forces  are  not  properly  evaluated  if  the  element  is 
skewed.  To  remedy  this  drawback,  Sa(0)  are  replaced  by  the  uniform  matrices,  Sa,  defined 
by  Belytschko  et  al.  [18,19]. 

^  I  BAi,ri,OdV  (159) 

where  Ve  is  the  element  volume.  Similarly,  equation  (85)  is  modified  as 


Ba  (0) 


bia 

b2a 

hsa 


and  the  stabilization  vectors  are  redefined  as 


(160) 


7a  —  (^h^Xi)  hi  (161) 

which  span  the  proper  null-space.  Since,  the  element  internal  force  vector  can  be  evaluated 
exactly  when  the  element  is  subjected  to  a  constant  strain  rate  field,  the  use  of  the  uniform 
gradient  matrices  Ba  leads  to  a  new  four-quadrature-point-element  which  passes  the  patch 
test. 
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4.  NUMERICAL  EXAMPLES 


We  now  present  geometrically  non-linear  analysis  of  flat  and  curved  composite  laminates 
and  compare  the  results  with  the  published  data  by  various  investigators.  This  task  is  per¬ 
formed  to  verify  and  validate  the  proposed  model  with  the  established  benchmark  problems. 
The  main  feature  of  the  proposed  formulation  is  that  it  can  easily  model  the  layered  material 
with  orthotropic  material  properties.  Comparison  is  made  with  some  of  the  celebrated  shell 
models  in  the  literature. 

4.1  Geometrically  nonlinear  bending  of  a  narrow  cantilever  plate  under  point 
load 

The  first  simulation  is  that  of  a  narrow  cantilever  isotropic  plate  subjected  to  a  concen¬ 
trated  load,  and  comparison  is  done  with  the  shell  element  SHEBA  by  Argyris  et  al.  [26]. 
The  concentrated  load  is  applied  in  equal  intervals  at  the  middle  of  the  right  boundary  (see 
Fig.  4).  The  values  of  Youngs  modulus  E  =  2.1x  10+®  kg/crri^  and  Poisson  ratio  v  =  0.3  are 
used  in  this  analysis.  The  value  of  applied  load,  P  is  4  x  10^  kg.  The  applied  load  versus  the 
displacement  path  (Fig.  5)  of  the  load  point  in  our  simulation  shows  a  close  agreement  with 
that  given  by  Argyris  et.  al.  [26].  The  superposed  initial  and  final  meshes  are  shown  in  Fig. 
6  and  surface  stresses  (Txx  are  shown  in  Fig.  7. 
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Displacement 


Nonlinear  bending  of  thin 
cantilever  isotropic  plate 


Fig.  4-  Schematic  diagram  of  the  problem. 

Nonlinear  bending  of  a 
thin  cantilever  isotropic  plate 


Fig.  5.  Load-deflection  diagram  at  point  under  the  load. 
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Fig.  7. 


Surface  stress  distribution  for  on  the  final  deformed  geometry. 
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4.2  Large  deformation  analysis  of  three-dimensional  curvilinear  beam 

This  example  presents  large  displacement  response  of  a  three-dimensional  curviUnear  can¬ 
tilever  beam.  The  beam  has  a  cross-sectional  area  of  1  in^  and  lies  in  the  X  —  Y  plane.  In  the 
undeformed  configuration  the  cantilever  bend  has  a  radius  of  100  inch  and  an  angle  of  45°. 
A  concentrated  load  is  applied  at  one  of  the  ends  of  the  beam  in  the  Z  direction  while  the 
opposite  end  is  kept  fixed.  Because  of  its  curved  geometry  and  the  eccentricity  of  the  applied 
load,  the  beam  undergoes  bending  as  well  as  twisting  deformations.  The  finite  element  mesh 
is  composed  of  16  x  2  x  1  elements.  Material  properties  used  are:  Young’s  modulus  E  =  lO'^psi 
and  Poisson’s  ratios  p  —  0. 

Figure  8  shows  the  initial  and  the  final  defiected  shapes  of  the  cantilever.  The  graph  for  tip 
deflection  plotted  against  the  load  parameter  can  be  seen  in  Figure  9.  Results  are  compared 
with  the  beam  element  of  Bathe  and  Bolourchi  and  3D  co-rotational  element  of  Moita  and 
Crisfield.  Computed  results  show  a  very  good  correlation  with  the  published  results. 


Nondimensionat  Displacement 


Fig.  8.  Initial  and  final  deformed  geometry. 
Three  Dimensional  Curviiinear  Beam 
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4.3  Geometrically  nonlinear  analysis  of  a  clamped  isotropic  plate  under  uniform 
pressure  load 

The  next  numerical  simulation  is  a  large  deformation  analysis  of  a  fully  clamped  isotropic 
plate.  Again  we  have  selected  this  isotropic  problem  to  benchmark  our  element  for  geomet¬ 
rically  nonlinear  analysis.  Material  properties  are,  Youngs  Modulus  E  =  2  x  10‘^kg/mm? 
and  Poisson  ratio  u  =  0.3.  The  plate  is  loaded  by  applying  a  uniform  pressure  of  intensity 
g  =  0.8  X  lO”**  kg/mm?  (see  Fig.  10).  To  solve  this  problem  we  generated  a  mesh  of  8  x  8  x  2 
hybrid  elements.  The  pseudo-time  increment  in  the  nonlinear  solution  strategy  is  At  =  0.01. 
A  total  of  100  steps  were  required  to  apply  the  total  load.  Fig  11  presents  the  nonlinear 
load-deflection  curve  obtained  at  the  point  under  the  load  which  compares  very  well  with  the 
experiment  load-deflection  curve  reported  in  Kawai  et  al.  [24].  Figs.  12  and  13  present  the 
top  and  bottom  view  of  the  surface  stresses  for  the  flnal  step. 
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Fig.  1 0.  Schematic  diagram  of  the  problem. 


Nonlinear  bending  of  isotropic  plate 
under  uniform  distributed  load 


Fig.  11.  Load-deflection  diagram  at  the  center  of  the  plate 
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STRESS  1 
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-1.07E+00 
-6.12E-01 
-1 .57E-01 
2.97E-01 
7.52E-01 
1.21E+00 
1.66E+00 

Current  View 
Min  =  -1.52E+00 
X  =  8.70E-03 

Y  =  2.00E+02 
Z  =  2.60E-01 
Max=  1.66E+00 
X  =  1 .84E-03 

Y  =  2.00E+02 
Z=-1.72E+00 


Fig.  12.  Surface  stress  distribution  for  Gxx  on  the  top  surface. 
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Current  View 
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Z=-1.72E+00 


Fig.  13.  Surface  stress  distribution  for  on  the  bottom  surface. 
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4.4  Pinched  cylindrical  shell  with  free  edges 


The  next  simulation  is  that  of  a  cylindrical  shell  with  free  edges  (Fig.  14),  subjected 
to  two  opposite  point  loads.  In  this  test  case  the  shell  is  able  to  undergo  finite  rotations 
and  thus  provides  a  severe  test  for  the  veracity  of  the  underlying  formulation.  The  overall 
response  represents  two  distinct  features;  an  initial  portion  that  is  dominated  by  bending  and 
is  characterized  by  large  displacement,  and  a  later  portion  that  is  dominated  by  membrane 
effects  that  are  characterized  by  a  very  stiff  response.  Invoking  the  symmetry  conditions,  only 
one-eighth  of  the  shell  is  modeled  using  16  x  8  x  2  elements  (16  along  the  periphery,  8  along 
the  length  and  2  through  the  thickness).  The  depicted  curves  show  the  displacement  under 
the  load,  as  well  as  the  horizontal  displacement  at  the  side  of  the  shell.  The  final  deformed 
configuration  of  the  shell  is  shown  in  Fig.  16. 
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Fig.  14-  Initial  geometry  and  loads. 
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point  A 
point  B 
#  Jiang  et.  al  point  A 

»  Jiang  et.  al  point  B 


Displacement  (m) 


Fig.  15.  Load-deflection  diagram. 
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Fig.  16.  Final  deformed  geometry. 


4.5  Pinched  hemispherical  shell 


The  next  problem  is  that  of  a  pinched  hemispherical  shell  which  is  a  popular  benchmark 
problem  for  hnear  as  well  as  nonlinear  shell  analysis.  The  undeformed  configuration  of  the 
shell  has  a  18°  hole  at  the  top  and  is  subjected  to  two  inward  and  two  outward  forces  that 
are  90°  apart.  The  material  and  geometric  properties  are:  E  =  6.825  x  psi,  u  =  0.3, 
radius  =  10  inch  and  thickness  t  =  0.04  inch.  Because  of  the  symmetry  conditions  that  can 
be  apphed,  only  one  quadrant  needs  to  be  modeled.  A  plot  of  the  pinching  load  versus  the 
deflection  under  the  corresponding  pinching  load  is  shown  in  Fig.  18.  Comparison  is  made 
with  “stress  resultant  geometrically  exact  shell  model  of  Simo  et  al.  25] .  Figure  19  shows  the 
final  deformed  mesh  configuration  without  any  magnification  of  the  deformation. 


Fig.  1 7.  Initial  geometry  and  loads. 
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Load 


•  Simoet.  al  x-dir 
^  Simoet.  al  y-dir 


Displacement 


Fig.  18.  Load-deflection  diagram. 
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Fig.  19.  Final  deformed  geometry. 


References 


[1]  Zienkiewicz,  O.C.,  (1977),  The  Finite  Element  Method  (3rd  edition),  McGraw-Hill. 

[2]  Bathe,  K.J.,  (1982),  Finite  Element  Procedures  in  Engineering  Analysis,  Prentice-Hall, 
Inc.  New  Jersey. 

[3]  Liao,  C.L.  and  Reddy,  J.N.,  (1990),  “Analysis  of  anisotropic,  stiffened  composite  laminates 
using  a  continuum-based  shell  element” ,  Comp.  Struct.,  34,  No.  6,  805-815. 

[4]  Belytschko,  T.  and  Hsieh,  B.J.,  (1973),  “Non-linear  transient  finite  element  analysis  with 
convected  coordinates”,  Intemat.  J.  Numer.  Methods  Engrg.,  7,  255-271. 

[5]  Wempner,  G.,  (1969),  “Finite  elements,  finite  rotations  and  small  strains  of  flexible  shells” , 
it  Internat.  J.  Solids  Structures,  5,  117-153. 

[6]  Harrigmoe,  G.  and  Bergan,  P.G.,  (1978),  “Nonhnear  analysis  of  free-form  shells  by  flat 
finite  elements”,  Comput.  Methods  in  Appl.  Mech.  Engrg.,  16,  11-35. 

[7]  Clarke,  M.J.  and  Hancock,  G.J.,  (1990),  “A  study  of  incremental-iterative  strategies  for 
non-linear  analysis”,  Intemat.  J.  Numer.  Methods  Engrg.,  29,  1365-1391. 

[8]  Bergan,  P.G.,  Harrigmoe,  G.,  Krakeland,  B.  and  Soreide,  T.H.,  (1978),  “Solution  tech¬ 
niques  for  non-linear  finite  element  problems”,  Intemat.  J.  Numer.  Methods  Engrg.,  12, 
1677-1696. 

[9]  Argyris,  J.,  Balmer,  H.  and  Doltsinis,  I.S.,  (1987),  “Implementation  of  a  nonlinear  capa¬ 
bility  on  a  hnear  software  system”,  Comput.  Methods  in  Appl.  Mech.  Engrg.,  65,  267-291. 

[10]  Riks,  E.,  (December  1971),  “The  application  of  Newton’s  method  to  the  problem  of  elastic 
stability”,  ASME,  J.  of  Appl.  Mech.,  1060-1065. 


3-49 


[11]  Riks,  E.,  (1979),  “An  incremental  approach  to  the  solution  of  snapping  and  buckling 
problems”,  Intemat.  J.  Solids  Structures,  15,  529-551. 

[12]  Crisfield,  M.A.,  (1981),  “A  fast  incremental/iterative  solution  procedure  that  handles 
’’snap-through””.  Comp.  Struct.,  13,  55-62. 

[13]  Bathe,  K.J.  and  Dvorkin,  E.N.,  (1983),  “On  the  automatic  solution  of  nonlinear  finite 
element  equations”,  Comp.  Struct.,  17,  No.  5-6,  871-879. 

[14]  Malkus,  D.S.  and  Hughes,  T.J.R.,  (1978),  “Mixed  finite  element  methods-reduced  and 
selective  integration  techniques:  a  unification  of  conceptes”,  Comput.  Methods  in  Appl. 
Mech.  Engrg.,  15,  63-81. 

[15]  Nagtegaal,  J.C.,  Parks,  D.M.  and  Rice,  J.R.,  (1974),  “On  numerical  accurate  finite  element 
solutions  in  the  fully  plastic  range”,  Comput.  Methods  in  Appl.  Mech.  Engrg.,  4,  153-177. 

[16]  Flanagan,  D.P.  and  Belytschko,  T.,  (1981),  “A  uniform  strain  hexahedron  and  quadri¬ 
lateral  with  orthogaonal  hourglass  control”,  Intemat.  J.  Numer.  Methods  Engrg.,  17, 
679-706. 

[17]  Belytschko,  T.,  (1983),  “Correction  of  article  by  Flanagan,  D.P.  and  Belytschko,  T.”, 
Intemat.  J.  Numer.  Methods  Engrg.,  19,  467-468. 

[18]  Belytschko,  T.,  Ong,  J.S.-J.,  Liu,  W.K.  and  Kennedy,  J.M.,  (1984),  “Hourglass  control  in 
linear  and  nonlinear  problems”,  Comput.  Methods  in  Appl.  Mech.  Engrg.,  43 , 251-276. 

[19]  Liu,  W.K.  and  Belytschko,  T.,  (1984),  “Efficient  linear  and  nonlinear  heat  condcution 
with  a  quadrilateral  element”,  Intemat.  J.  Numer.  Methods  Engrg.,  20,  931-948. 

[20]  Belytshcko,  T.  and  Bindeman,  L.P.,  (1991),  “Assumed  strain  stabifization  of  the  4-node 
quadrilateral  with  1-point  quadrature  for  nonlinear  problems” ,  Comput.  Methods  in  Appl. 
Mech.  Engrg.,  88,  311-340. 


3-50 


[21]  Belytshcko,  T.  and  Bindeman,  L.P.,  “Assumed  strain  stabilization  of  the  eight  node  hex- 
ahedral  element”,  to  be  published. 

[22]  Liu,  W.K.,  Ong,  J.S.-J.  and  Uras,  R.A.,  (1985),  “Finite  element  stabilization  matriced  - 
a  unification  approach”,  Comput.  Methods  in  Appl.  Mech.  Engrg.,  53,  13-46. 

[23]  Liu,  W.K.,  Chang,  H.,  Chen,  J.S.  and  Belytschko,  T.,  (1988),  “ALE  Petrov- Galerkin  finite 
elements  for  nonlinear. continua”,  CoTnput.  Methods  in  Appl.  Mech.  Engvg..,  68,  259-310. 

[24]  Kawai,  T.  and  Yoshimura,  N.,  (1969)  “Analysis  of  large  deflection  of  plates  by  the  finite 
element  rnethod”,  Intemat.  J.  Numer.  Methods  Engrg.  1,  123-133. 

[25]  Simo,  J.C.,  Fox,  D.D  and  Rifai,  M.S.,  (1990)  “On  a  stress  resultant  geometrically  exact 
shell  model.  Part  III:  Computational  aspects  of  the  nonlinear  theory,”  Comput.  Methods 
in  Appl.  Mech.  Engrg.,  79  ,  21-70. 

[26]  Argyris,  J.  and  Tenek,  L.,  (1994)  “An  efficient  and  locking-free  flat  anisotropic  plate  and 
shell  triangular  element”,  Comput  Methods  in  Appl.  Mech.  Engrg.,  118,  63-119. 

[27]  Jiang,  L.,  Chemuka,  M.W.,  and  Pegg,  N.G.  ,  (1994)  “  A  co-rotational  updated  Lagrangiait 
formulation  for  geometrically  nonlinear  finite  element  analysis  of  shell  structures  ,  Finite 
Element  in  Analysis  and  design,  Vol.  18, 129-140. 


3-51 


Part  IV 


Simulation  of  Coupled  Systems 


Chapter  1 

Coupled  Elastic/SMA  Simulations 


This  section  is  dedicated  to  numerical  simulations  of  coupled  elastic-shape  memory  alloy 
composite  systems.  The  simulations  depict  an  initial  configuration  is  given  which  is  then 
transformed  via  a  thermal  variation  to  a  final  configuration.  Several  simulations  are  presented 
to  shown  the  robustness  of  the  shape  memory  alloy  and  its  capacity  to  be  embedded  within 
an  elastic  medium  to  enable  a  global  shape  changes. 
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1.1  Two-Dimensional  Cantilever  Beam 


The  first  simulation  consists  of  a  configuration  with  an  initially  curved  two  dimensional 
cantilever  beam.  The  initial  configuration  of  the  beam  is  shown  in  Figure  1.1.  The  simulation 
is  conducted  with  a  total  net  reaction  of  the  beam  of  zero  and  an  initial  thermal  state  of 
31  °C  which  is  then  heated  to  54  "C.  The  material  under  consideration  is  NiTi  and  its 
associated  material  parameters  can  be  found  in  Table  1.1.  The  progression  of  the  simulation 
is  shown  in  Figures  1.2  and  1.3. 
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Figure  1.1:  Initial  Configuration  of  the  Composite  Elastic/SMA  Beam. 


Table  1.1:  Material  Properties  for  NiTi.(BRiNsoN  Lammering  [1993]) 
Young’s  Moduli  Em  =  Ea  =  67  GPa 

Critical  stresses  for  de- twinning  =  100  MPa  and  =  170  MPa 

Martensite  production  temperatures  =  18.4  C  and  Tmf  =  9  C 
Austenite  production  temperatures  To,  =  34.5  C  and  Taf  =  49  C 
Austenite  production  slope  Ca  =  13.8  MPa/C 
Martensite  production  slope  Cm  =  8  MPa/C 
Maximum  transformation  strain  =  0.067 

Thermal  expansion  coefficient  a  =  6.5  /istrain/C _ 
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T  =  T  =  5Z.75°C 


Figure  1.2;  Simulation  of  an  elastic/sma  beam  under  a  thermal  cycle. 


1.2  Two-Dimensional  Simply  Supported  Beam 


The  second  simulation  consists  of  a  configuration  with  an  initially  curved  two  dimensional 
simply  supported  beam.  The  initial  configuration  of  the  beam  is  shown  in  Figure  1.4.  The 
simulation  is  conducted  with  a  total  net  reaction  of  the  beam  of  zero  and  an  initial  thermal 
state  of  36.3  °C  which  is  then  heated  to  50.3  °C.  The  material  under  consideration  is  NiTi 
and  its  associated  material  parameters  can  be  found  in  Table  1.2.  The  progression  of  the 
simulation  is  shown  in  Figures  1.5  and  1.6. 
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Figure  1.4:  Initial  Configuration  of  the  Composite  Elastic/SMA  Beam. 


Table  1.2:  Material  Properties  for  NiTi. (Brinson  fc  Lammering  [1993]) 
Young’s  Moduli  =  Ea  =  67  GPa 

Critical  stresses  for  de- twinning  <7*^  =  100  MPa  and  =  170  MPa 
Martensite  production  temperatures  Tma  =  18.4  C  and  T^f  =  9  C 
Austenite  production  temperatures  Taa  =  34.5  C  and  Taf  =  49  C 
Austenite  production  slope  Ca  =  13.8  MPa/C 
Martensite  production  slope  Cm  =  8  MPa/C 
Maximum  transformation  strain  ei  =  0.067 
Thermal  expansion  coefficient  a  =  6.5  ^strain/C 
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